library(tidyverse) # for data wrangling etc
library(rstanarm) # for fitting models in STAN
library(cmdstanr) # for cmdstan
library(brms) # for fitting models in STAN
library(standist) # for exploring distributions
library(HDInterval) # for HPD intervals
library(posterior) # for posterior draws
library(coda) # for diagnostics
library(bayesplot) # for diagnostics
library(ggmcmc) # for MCMC diagnostics
library(rstan) # for interfacing with STAN
library(emmeans) # for marginal means etc
library(broom) # for tidying outputs
library(DHARMa) # for residual diagnostics
library(tidybayes) # for more tidying outputs
library(ggeffects) # for partial plots
library(broom.mixed) # for tidying MCMC outputs
library(patchwork) # for multiple plots
library(ggridges) # for ridge plots
library(bayestestR) # for ROPE
library(see) # for some plots
library(easystats) # framework for stats, modelling and visualisation
library(modelsummary)
source("helperFunctions.R")Bayesian GLMM Part2
1 Preparations
Load the necessary libraries
2 Scenario
To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.
A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.
| FISHID | MASS | TRIAL | SMR_contr | CHANGE |
|---|---|---|---|---|
| 1 | 35.69 | LowSalinity | 5.85 | -31.92 |
| 2 | 33.84 | LowSalinity | 6.53 | 2.52 |
| 3 | 37.78 | LowSalinity | 5.66 | -6.28 |
| .. | .. | .. | .. | .. |
| 1 | 36.80 | HighTemperature | 5.85 | 18.32 |
| 2 | 34.98 | HighTemperature | 6.53 | 19.06 |
| 3 | 38.38 | HighTemperature | 5.66 | 19.03 |
| .. | .. | .. | .. | .. |
| 1 | 45.06 | Hypoxia | 5.85 | -18.61 |
| 2 | 43.51 | Hypoxia | 6.53 | -5.37 |
| 3 | 45.11 | Hypoxia | 5.66 | -13.95 |
| FISHID | Categorical listing of the individual fish that are repeatedly sampled |
| MASS | Mass (g) of barramundi. Covariate in analysis |
| TRIAL | Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen. |
| SMR_contr | Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia) |
| CHANGE | Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'. |
3 Read in the data
Rows: 180 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): TRIAL
dbl (4): FISHID, MASS, SMR_contr, CHANGE
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 180
Columns: 5
$ FISHID <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
$ MASS <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
$ TRIAL <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
$ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
$ CHANGE <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…
# A tibble: 6 × 5
FISHID MASS TRIAL SMR_contr CHANGE
<dbl> <dbl> <chr> <dbl> <dbl>
1 1 35.7 LowSalinity 5.85 -31.9
2 2 33.8 LowSalinity 6.53 2.52
3 3 37.8 LowSalinity 5.66 -6.28
4 4 26.6 LowSalinity 6.28 -4.35
5 5 37.6 LowSalinity 4.41 -3.07
6 6 37.7 LowSalinity 4.82 -15.1
spc_tbl_ [180 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ FISHID : num [1:180] 1 2 3 4 5 6 7 8 9 10 ...
$ MASS : num [1:180] 35.7 33.8 37.8 26.6 37.6 ...
$ TRIAL : chr [1:180] "LowSalinity" "LowSalinity" "LowSalinity" "LowSalinity" ...
$ SMR_contr: num [1:180] 5.85 6.53 5.66 6.28 4.41 ...
$ CHANGE : num [1:180] -31.92 2.52 -6.28 -4.35 -3.07 ...
- attr(*, "spec")=
.. cols(
.. FISHID = col_double(),
.. MASS = col_double(),
.. TRIAL = col_character(),
.. SMR_contr = col_double(),
.. CHANGE = col_double()
.. )
- attr(*, "problems")=<externalptr>
norin (180 rows and 5 variables, 5 shown)
ID | Name | Type | Missings | Values | N
---+-----------+-----------+----------+------------------+-----------
1 | FISHID | numeric | 0 (0.0%) | [1, 60] | 180
---+-----------+-----------+----------+------------------+-----------
2 | MASS | numeric | 0 (0.0%) | [22.52, 64.07] | 180
---+-----------+-----------+----------+------------------+-----------
3 | TRIAL | character | 0 (0.0%) | HighTemperature | 60 (33.3%)
| | | | Hypoxia | 60 (33.3%)
| | | | LowSalinity | 60 (33.3%)
---+-----------+-----------+----------+------------------+-----------
4 | SMR_contr | numeric | 0 (0.0%) | [3.98, 6.7] | 180
---+-----------+-----------+----------+------------------+-----------
5 | CHANGE | numeric | 0 (0.0%) | [-55.01, 116.33] | 180
---------------------------------------------------------------------
| Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|
| FISHID | 60 | 0 | 30.5 | 17.4 | 1.0 | 30.5 | 60.0 | |
| MASS | 175 | 0 | 40.2 | 8.0 | 22.5 | 39.3 | 64.1 | |
| SMR_contr | 60 | 0 | 5.2 | 0.6 | 4.0 | 5.1 | 6.7 | |
| CHANGE | 180 | 0 | 19.6 | 33.8 | -55.0 | 16.4 | 116.3 | |
| TRIAL | N | % | ||||||
| HighTemperature | 60 | 33.3 | ||||||
| Hypoxia | 60 | 33.3 | ||||||
| LowSalinity | 60 | 33.3 |
| TRIAL | Unique | Missing Pct. | Mean | SD | Min | Median | Max | Histogram | |
|---|---|---|---|---|---|---|---|---|---|
| FISHID | LowSalinity | 60 | 0 | 30.5 | 17.5 | 1.0 | 30.5 | 60.0 | |
| HighTemperature | 60 | 0 | 30.5 | 17.5 | 1.0 | 30.5 | 60.0 | ||
| Hypoxia | 60 | 0 | 30.5 | 17.5 | 1.0 | 30.5 | 60.0 | ||
| MASS | LowSalinity | 60 | 0 | 38.1 | 7.2 | 23.8 | 37.6 | 53.0 | |
| HighTemperature | 60 | 0 | 38.3 | 7.2 | 22.5 | 37.7 | 53.6 | ||
| Hypoxia | 60 | 0 | 44.2 | 8.1 | 28.5 | 44.0 | 64.1 | ||
| SMR_contr | LowSalinity | 60 | 0 | 5.2 | 0.6 | 4.0 | 5.1 | 6.7 | |
| HighTemperature | 60 | 0 | 5.2 | 0.6 | 4.0 | 5.1 | 6.7 | ||
| Hypoxia | 60 | 0 | 5.2 | 0.6 | 4.0 | 5.1 | 6.7 | ||
| CHANGE | LowSalinity | 60 | 0 | 8.2 | 34.4 | -55.0 | 6.4 | 99.3 | |
| HighTemperature | 60 | 0 | 50.4 | 22.8 | 7.0 | 52.2 | 116.3 | ||
| Hypoxia | 60 | 0 | 0.2 | 17.1 | -41.2 | -2.7 | 46.6 | ||
| TRIAL | N | % | |||||||
| HighTemperature | 60 | 33.3 | |||||||
| Hypoxia | 60 | 33.3 | |||||||
| LowSalinity | 60 | 33.3 |
4 Exploratory data analysis
Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 \sim{} \mathcal{N}(16, 35)\\ \beta_{1-6} \sim{} \mathcal{N}(0, 70)\\ \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.
Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.
Conclusions:
- these seem reasonable enough
The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.
Doing so would introduce two additional assumptions:
- linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
- homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth(method = "lm") +
geom_point()`geom_smooth()` using formula = 'y ~ x'
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth() +
geom_point()`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.
Conclusions:
- generally, the difference between control and trial SMR is greatest for the High temperature trial.
- there is less consistency about the relative effects of Hypoxia and Low Salinity.
- hence a random intercept/slope model might be useful.
Finally, as an acknowledgement that the metabolic response might be influenced
by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.
- add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardised by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
- add them as a regular covariate. This will include the covariate in the model and estimate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.
Lets explore the relationship between the response and Mass separately for each Trial.
`geom_smooth()` using formula = 'y ~ x'
Conclusions:
- there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
- there is no evidence that the response and Mass are related in a non-linear manner.
5 Fit the model
In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.
Priors for model 'norin.rstanarm'
------
Intercept (after predictors centered)
Specified prior:
~ normal(location = 20, scale = 2.5)
Adjusted prior:
~ normal(location = 20, scale = 85)
Coefficients
Specified prior:
~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
Adjusted prior:
~ normal(location = [0,0,0,...], scale = [178.99,178.99,135.15,...])
Auxiliary (sigma)
Specified prior:
~ exponential(rate = 1)
Adjusted prior:
~ exponential(rate = 0.03)
Covariance
~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more details
This tells us:
- for the intercept (when the family is Gaussian), a normal prior with a mean of 20 and a standard deviation of 85 is used. A value of 2.5 is used for scaling all parameter standard deviations. The value of 20 is based on the mean of the response (
mean(norin$CHANGE)) and the scaled standard deviation of 85 is based on multiplying the scaling factor by the standard deviation of the response.
- for the coefficients (in this case, just the difference between strong and weak innoculation), the default prior is a normal prior centered around 0 with a standard deviations of 178.99, 178.99, 135.15, 10.60, 34.2 and 34.2 repectively. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the stanard devation of the numerical dummy variables for the predictor (then rounded).
(Intercept) TRIALHypoxia
Inf 178.99305
TRIALLowSalinity SMR_contr
178.99305 135.14517
MASS TRIALHypoxia:SMR_contr
10.59372 34.19156
TRIALLowSalinity:SMR_contr
34.19156
- the auxillary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be \(\sigma\), for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in
rstanarm, this is a exponential with a rate of 1 which is then adjusted by devision with the standard deviation of the response.
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Conclusions:
- we see that the range of predictions is faily wide and the predicted means could range from a small negative number to a relatively large positive number.
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
- \(\beta_0\): normal centered at 17 with a standard deviation of 35
- mean of 17: since
median(norin$CHANGE) - sd of 35: since
mad(norin$CHANGE)
- mean of 17: since
- \(\beta_1\): normal centred at 0 with a standard deviation of 70
- sd of 70: since
sd(norin$CHANGE) / apply(model.matrix(~scale(SMR_contr)*TRIAL + scale(MASS), norin), 2, sd)
- sd of 70: since
- \(\sigma\): exponential with rate of 0.03
- since
1 / sd(norin$CHANGE)
- since
- \(\Sigma\): decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.
I will also overlay the raw data for comparison.
norin.rstanarm2 <- stan_glmer(CHANGE ~ (1 | FISHID) + TRIAL * scale(SMR_contr, scale = FALSE) + offset(MASS),
data = norin,
family = gaussian(),
prior_intercept = normal(17, 35, autoscale = FALSE),
prior = normal(0, 70, autoscale = FALSE),
prior_aux = rstanarm::exponential(0.03, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
prior_PD = TRUE,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Now lets refit, conditioning on the data.
norin.rstanarm4 <- stan_glmer(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + offset(MASS),
data = norin,
family = gaussian(),
prior_intercept = normal(17, 35, autoscale = FALSE),
prior = normal(0, 70, autoscale = FALSE),
prior_aux = rstanarm::exponential(0.03, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
preds <- norin.rstanarm4 |> posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE
)
plot(norin.resids, quantreg = FALSE) [1] "(Intercept)"
[2] "TRIALHypoxia"
[3] "TRIALLowSalinity"
[4] "scale(SMR_contr, scale = FALSE)"
[5] "TRIALHypoxia:scale(SMR_contr, scale = FALSE)"
[6] "TRIALLowSalinity:scale(SMR_contr, scale = FALSE)"
[7] "b[(Intercept) FISHID:1]"
[8] "b[(Intercept) FISHID:2]"
[9] "b[(Intercept) FISHID:3]"
[10] "b[(Intercept) FISHID:4]"
[11] "b[(Intercept) FISHID:5]"
[12] "b[(Intercept) FISHID:6]"
[13] "b[(Intercept) FISHID:7]"
[14] "b[(Intercept) FISHID:8]"
[15] "b[(Intercept) FISHID:9]"
[16] "b[(Intercept) FISHID:10]"
[17] "b[(Intercept) FISHID:11]"
[18] "b[(Intercept) FISHID:12]"
[19] "b[(Intercept) FISHID:13]"
[20] "b[(Intercept) FISHID:14]"
[21] "b[(Intercept) FISHID:15]"
[22] "b[(Intercept) FISHID:16]"
[23] "b[(Intercept) FISHID:17]"
[24] "b[(Intercept) FISHID:18]"
[25] "b[(Intercept) FISHID:19]"
[26] "b[(Intercept) FISHID:20]"
[27] "b[(Intercept) FISHID:21]"
[28] "b[(Intercept) FISHID:22]"
[29] "b[(Intercept) FISHID:23]"
[30] "b[(Intercept) FISHID:24]"
[31] "b[(Intercept) FISHID:25]"
[32] "b[(Intercept) FISHID:26]"
[33] "b[(Intercept) FISHID:27]"
[34] "b[(Intercept) FISHID:28]"
[35] "b[(Intercept) FISHID:29]"
[36] "b[(Intercept) FISHID:30]"
[37] "b[(Intercept) FISHID:31]"
[38] "b[(Intercept) FISHID:32]"
[39] "b[(Intercept) FISHID:33]"
[40] "b[(Intercept) FISHID:34]"
[41] "b[(Intercept) FISHID:35]"
[42] "b[(Intercept) FISHID:36]"
[43] "b[(Intercept) FISHID:37]"
[44] "b[(Intercept) FISHID:38]"
[45] "b[(Intercept) FISHID:39]"
[46] "b[(Intercept) FISHID:40]"
[47] "b[(Intercept) FISHID:41]"
[48] "b[(Intercept) FISHID:42]"
[49] "b[(Intercept) FISHID:43]"
[50] "b[(Intercept) FISHID:44]"
[51] "b[(Intercept) FISHID:45]"
[52] "b[(Intercept) FISHID:46]"
[53] "b[(Intercept) FISHID:47]"
[54] "b[(Intercept) FISHID:48]"
[55] "b[(Intercept) FISHID:49]"
[56] "b[(Intercept) FISHID:50]"
[57] "b[(Intercept) FISHID:51]"
[58] "b[(Intercept) FISHID:52]"
[59] "b[(Intercept) FISHID:53]"
[60] "b[(Intercept) FISHID:54]"
[61] "b[(Intercept) FISHID:55]"
[62] "b[(Intercept) FISHID:56]"
[63] "b[(Intercept) FISHID:57]"
[64] "b[(Intercept) FISHID:58]"
[65] "b[(Intercept) FISHID:59]"
[66] "b[(Intercept) FISHID:60]"
[67] "sigma"
[68] "Sigma[FISHID:(Intercept),(Intercept)]"
[69] "accept_stat__"
[70] "stepsize__"
[71] "treedepth__"
[72] "n_leapfrog__"
[73] "divergent__"
[74] "energy__"
posterior_vs_prior(norin.rstanarm3,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y"),
regex_pars = "^.Intercept|TRIAL|SMR_contr|MASS|sigma|Sigma"
)
Drawing from prior...
Conclusions:
- in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.
Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
'dim(x)[MARGIN]'
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.
Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.
norin.form <- bf(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+offset(MASS),
family = gaussian()
)
options(width=100)
norin.form |> get_prior(data=norin) prior class coef group resp dpar nlpar
(flat) b
(flat) b SMR_contr
(flat) b TRIALHypoxia
(flat) b TRIALHypoxia:SMR_contr
(flat) b TRIALLowSalinity
(flat) b TRIALLowSalinity:SMR_contr
student_t(3, -23.7797222222222, 34.5) Intercept
student_t(3, 0, 34.5) sd
student_t(3, 0, 34.5) sd FISHID
student_t(3, 0, 34.5) sd Intercept FISHID
student_t(3, 0, 34.5) sigma
lb ub source
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
0 default
0 (vectorized)
0 (vectorized)
0 default
norin.brm <- brm(norin.form,
data=norin,
iter = 5000,
warmup = 1000,
chains = 3, cores = 3,
thin = 5,
refresh = 0,
backend = "rstan")Compiling Stan program...
Start sampling
This tells us:
- for the intercept (when the family is Gaussian), a t-distribution prior with a mean of 16.4 and a standard deviation of 34.5 is used:
- mean of 16.4, since
median(norin$CHANGE) - sd of 34.5 since
mad(norin$CHANGE)
- mean of 16.4, since
brmsuses flat (inproper) priors for all population effects- the prior on sigma is a half t-distribution with a mean of 0 and standard deviation of 34.5:
- sd of 34.5 since
mad(norin$CHANGE)
- sd of 34.5 since
- the hyperprior on the standard deviation is a half t-distribution with a mean of 0 and a standard deviation of 34.5.
The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]
When defining our own priors, we typically do not want them to be scaled.
If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):
- \(\beta_0\): normal centred at 16 with a standard deviation of 35
- mean of 16: since
median(norin$CHANGE) - sd of 35: since
mad(norin$CHANGE)
- mean of 16: since
- \(\beta_{1-2}\) (TRIAL): normal centred at 0 with a standard deviation of 70
- sd of 70: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 70: since
- \(\beta_{3}\) (SMR_contr): normal centred at 0 with a standard deviation of 54
- sd of 54: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 54: since
- \(\beta_{4}\) (MASS): normal centred at 0 with a standard deviation of 4
- sd of 4: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 4: since
- \(\beta_{5-6}\) (interaction): normal centred at 0 with a standard deviation of 14
- sd of 14: since
sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
- sd of 14: since
- \(\sigma\): half-cauchy with parameters 0 and 6 or else a student_t with parameters 0 and 34
- since
sd(norin$CHANGE)
- since
- \(\sigma_j\): half-cauchy with parameters 0 and 5.8
- since
sqrt(sd(norin$CHANGE)) - we want this prior to have most mass close to zero for the purpose of regularisation
- since
- \(\Sigma\): decov with:
- regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
- concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
- shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.
Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:
- half-t: as the degrees of freedom approach infinity, this will approach a half-normal
- half-cauchy: this is essentially a half-t with 1 degree of freedom
- exponential
I will also overlay the raw data for comparison.
# A tibble: 3 × 3
TRIAL `median(CHANGE)` `mad(CHANGE)`
<fct> <dbl> <dbl>
1 HighTemperature 52.2 23.2
2 Hypoxia -2.66 13.7
3 LowSalinity 6.43 38.4
norin |>
group_by(TRIAL, FISHID) |>
summarise(
median = median(CHANGE),
MAD = mad(CHANGE)
) |>
ungroup(FISHID) |>
summarise(sd(median))`summarise()` has grouped output by 'TRIAL'. You can override using the
`.groups` argument.
# A tibble: 3 × 2
TRIAL `sd(median)`
<fct> <dbl>
1 HighTemperature 22.8
2 Hypoxia 17.1
3 LowSalinity 34.4
TRIALHypoxia TRIALLowSalinity
71.597218 71.597218
SMR_contr MASS
54.058070 4.237486
TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr
13.676626 13.676626
standist::visualize("gamma(2, 1)", "gamma(35, 1)",
"student_t(3,0, 34)",
"cauchy(0, 5.8)",
xlim = c(-10, 100)
)priors <- prior(normal(16, 35), class = "Intercept") +
prior(normal(0, 70), class = "b", coef = "TRIALHypoxia") +
prior(normal(0, 70), class = "b", coef = "TRIALLowSalinity") +
prior(normal(0, 54), class = "b", coef = "SMR_contr") +
prior(normal(0, 4), class = "b", coef = "MASS") +
prior(normal(0, 13), class = "b") +
prior(student_t(3, 0, 34), class = "sigma") +
## prior(gamma(35, 1), class = 'sigma') +
prior(cauchy(0, 5.8), class = "sd")
norin.form <- bf(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
family = gaussian()
)
norin.brm2 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3, cores = 3,
thin = 5,
refresh = 0,
backend = "rstan"
)Compiling Stan program...
Start sampling
Note in the above model, the output may have included a warning message alerting us the presence of divergent transitions. Divergent transitions are an indication that the sampler has encountered poor sampling conditions - the more divergent transitions, the more severe the issue.
Typically, divergent transitions are the result of either:
- a miss-specified model
- priors that permit the sampler to drift into unsupported areas
- complex posterior “features” (with high degrees of curvature) for which the sampler was inadequately tuned during the warmup phase
Accordingly, these divergent transitions can be addressed by either:
- reviewing the model structure
- adopting tighter priors
- increase the adaptive delta from the default of 0.8 to closer to 1. The adaptive delta defines the average acceptance probability that the sampler should aspire to during the warmup phase. Increasing the adaptive delta results in a smaller step size (and thus fewer divergences and more robust samples) however it will also result in slower sampling speeds.
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
norin.brm3 <- update(norin.brm2,
sample_prior = "yes",
control = list(adapt_delta = 0.99),
refresh = 0,
cores = 3
)The desired updates require recompiling the model
Compiling Stan program...
Start sampling
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
[1] "b_Intercept" "b_TRIALHypoxia"
[3] "b_TRIALLowSalinity" "b_SMR_contr"
[5] "b_MASS" "b_TRIALHypoxia:SMR_contr"
[7] "b_TRIALLowSalinity:SMR_contr" "sd_FISHID__Intercept"
[9] "sigma" "Intercept"
[11] "r_FISHID[1,Intercept]" "r_FISHID[2,Intercept]"
[13] "r_FISHID[3,Intercept]" "r_FISHID[4,Intercept]"
[15] "r_FISHID[5,Intercept]" "r_FISHID[6,Intercept]"
[17] "r_FISHID[7,Intercept]" "r_FISHID[8,Intercept]"
[19] "r_FISHID[9,Intercept]" "r_FISHID[10,Intercept]"
[21] "r_FISHID[11,Intercept]" "r_FISHID[12,Intercept]"
[23] "r_FISHID[13,Intercept]" "r_FISHID[14,Intercept]"
[25] "r_FISHID[15,Intercept]" "r_FISHID[16,Intercept]"
[27] "r_FISHID[17,Intercept]" "r_FISHID[18,Intercept]"
[29] "r_FISHID[19,Intercept]" "r_FISHID[20,Intercept]"
[31] "r_FISHID[21,Intercept]" "r_FISHID[22,Intercept]"
[33] "r_FISHID[23,Intercept]" "r_FISHID[24,Intercept]"
[35] "r_FISHID[25,Intercept]" "r_FISHID[26,Intercept]"
[37] "r_FISHID[27,Intercept]" "r_FISHID[28,Intercept]"
[39] "r_FISHID[29,Intercept]" "r_FISHID[30,Intercept]"
[41] "r_FISHID[31,Intercept]" "r_FISHID[32,Intercept]"
[43] "r_FISHID[33,Intercept]" "r_FISHID[34,Intercept]"
[45] "r_FISHID[35,Intercept]" "r_FISHID[36,Intercept]"
[47] "r_FISHID[37,Intercept]" "r_FISHID[38,Intercept]"
[49] "r_FISHID[39,Intercept]" "r_FISHID[40,Intercept]"
[51] "r_FISHID[41,Intercept]" "r_FISHID[42,Intercept]"
[53] "r_FISHID[43,Intercept]" "r_FISHID[44,Intercept]"
[55] "r_FISHID[45,Intercept]" "r_FISHID[46,Intercept]"
[57] "r_FISHID[47,Intercept]" "r_FISHID[48,Intercept]"
[59] "r_FISHID[49,Intercept]" "r_FISHID[50,Intercept]"
[61] "r_FISHID[51,Intercept]" "r_FISHID[52,Intercept]"
[63] "r_FISHID[53,Intercept]" "r_FISHID[54,Intercept]"
[65] "r_FISHID[55,Intercept]" "r_FISHID[56,Intercept]"
[67] "r_FISHID[57,Intercept]" "r_FISHID[58,Intercept]"
[69] "r_FISHID[59,Intercept]" "r_FISHID[60,Intercept]"
[71] "prior_Intercept" "prior_b_TRIALHypoxia"
[73] "prior_b_TRIALLowSalinity" "prior_b_SMR_contr"
[75] "prior_b_MASS" "prior_b_TRIALHypoxia:SMR_contr"
[77] "prior_b_TRIALLowSalinity:SMR_contr" "prior_sigma"
[79] "prior_sd_FISHID" "lprior"
[81] "lp__" "accept_stat__"
[83] "stepsize__" "treedepth__"
[85] "n_leapfrog__" "divergent__"
[87] "energy__"
While we are here, we might like to explore a random intercept/slope model
priors <- prior(normal(16, 35), class = "Intercept") +
prior(normal(0, 70), class = "b", coef = "TRIALHypoxia") +
prior(normal(0, 70), class = "b", coef = "TRIALLowSalinity") +
prior(normal(0, 54), class = "b", coef = "SMR_contr") +
prior(normal(0, 4), class = "b", coef = "MASS") +
## prior(normal(0, 70), class = 'b') +
prior(student_t(3, 0, 34), class = "sigma") +
prior(cauchy(0, 5.8), class = "sd") +
prior(lkj_corr_cholesky(1), class = "cor")
norin.form <- bf(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS,
## sigma ~ TRIAL*SMR_contr + offset(MASS) + (1|FISHID),
family = gaussian()
)
norin.brm4 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 1000,
chains = 3, cores = 3,
thin = 5,
refresh = 0,
control = list(adapt_delta = 0.99),
backend = "rstan"
)Compiling Stan program...
Start sampling
Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low
Warning: Examine the pairs() plot to diagnose sampling problems
[1] "b_Intercept"
[2] "b_TRIALHypoxia"
[3] "b_TRIALLowSalinity"
[4] "b_SMR_contr"
[5] "b_MASS"
[6] "b_TRIALHypoxia:SMR_contr"
[7] "b_TRIALLowSalinity:SMR_contr"
[8] "sd_FISHID__Intercept"
[9] "sd_FISHID__TRIALHypoxia"
[10] "sd_FISHID__TRIALLowSalinity"
[11] "cor_FISHID__Intercept__TRIALHypoxia"
[12] "cor_FISHID__Intercept__TRIALLowSalinity"
[13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
[14] "sigma"
[15] "Intercept"
[16] "r_FISHID[1,Intercept]"
[17] "r_FISHID[2,Intercept]"
[18] "r_FISHID[3,Intercept]"
[19] "r_FISHID[4,Intercept]"
[20] "r_FISHID[5,Intercept]"
[21] "r_FISHID[6,Intercept]"
[22] "r_FISHID[7,Intercept]"
[23] "r_FISHID[8,Intercept]"
[24] "r_FISHID[9,Intercept]"
[25] "r_FISHID[10,Intercept]"
[26] "r_FISHID[11,Intercept]"
[27] "r_FISHID[12,Intercept]"
[28] "r_FISHID[13,Intercept]"
[29] "r_FISHID[14,Intercept]"
[30] "r_FISHID[15,Intercept]"
[31] "r_FISHID[16,Intercept]"
[32] "r_FISHID[17,Intercept]"
[33] "r_FISHID[18,Intercept]"
[34] "r_FISHID[19,Intercept]"
[35] "r_FISHID[20,Intercept]"
[36] "r_FISHID[21,Intercept]"
[37] "r_FISHID[22,Intercept]"
[38] "r_FISHID[23,Intercept]"
[39] "r_FISHID[24,Intercept]"
[40] "r_FISHID[25,Intercept]"
[41] "r_FISHID[26,Intercept]"
[42] "r_FISHID[27,Intercept]"
[43] "r_FISHID[28,Intercept]"
[44] "r_FISHID[29,Intercept]"
[45] "r_FISHID[30,Intercept]"
[46] "r_FISHID[31,Intercept]"
[47] "r_FISHID[32,Intercept]"
[48] "r_FISHID[33,Intercept]"
[49] "r_FISHID[34,Intercept]"
[50] "r_FISHID[35,Intercept]"
[51] "r_FISHID[36,Intercept]"
[52] "r_FISHID[37,Intercept]"
[53] "r_FISHID[38,Intercept]"
[54] "r_FISHID[39,Intercept]"
[55] "r_FISHID[40,Intercept]"
[56] "r_FISHID[41,Intercept]"
[57] "r_FISHID[42,Intercept]"
[58] "r_FISHID[43,Intercept]"
[59] "r_FISHID[44,Intercept]"
[60] "r_FISHID[45,Intercept]"
[61] "r_FISHID[46,Intercept]"
[62] "r_FISHID[47,Intercept]"
[63] "r_FISHID[48,Intercept]"
[64] "r_FISHID[49,Intercept]"
[65] "r_FISHID[50,Intercept]"
[66] "r_FISHID[51,Intercept]"
[67] "r_FISHID[52,Intercept]"
[68] "r_FISHID[53,Intercept]"
[69] "r_FISHID[54,Intercept]"
[70] "r_FISHID[55,Intercept]"
[71] "r_FISHID[56,Intercept]"
[72] "r_FISHID[57,Intercept]"
[73] "r_FISHID[58,Intercept]"
[74] "r_FISHID[59,Intercept]"
[75] "r_FISHID[60,Intercept]"
[76] "r_FISHID[1,TRIALHypoxia]"
[77] "r_FISHID[2,TRIALHypoxia]"
[78] "r_FISHID[3,TRIALHypoxia]"
[79] "r_FISHID[4,TRIALHypoxia]"
[80] "r_FISHID[5,TRIALHypoxia]"
[81] "r_FISHID[6,TRIALHypoxia]"
[82] "r_FISHID[7,TRIALHypoxia]"
[83] "r_FISHID[8,TRIALHypoxia]"
[84] "r_FISHID[9,TRIALHypoxia]"
[85] "r_FISHID[10,TRIALHypoxia]"
[86] "r_FISHID[11,TRIALHypoxia]"
[87] "r_FISHID[12,TRIALHypoxia]"
[88] "r_FISHID[13,TRIALHypoxia]"
[89] "r_FISHID[14,TRIALHypoxia]"
[90] "r_FISHID[15,TRIALHypoxia]"
[91] "r_FISHID[16,TRIALHypoxia]"
[92] "r_FISHID[17,TRIALHypoxia]"
[93] "r_FISHID[18,TRIALHypoxia]"
[94] "r_FISHID[19,TRIALHypoxia]"
[95] "r_FISHID[20,TRIALHypoxia]"
[96] "r_FISHID[21,TRIALHypoxia]"
[97] "r_FISHID[22,TRIALHypoxia]"
[98] "r_FISHID[23,TRIALHypoxia]"
[99] "r_FISHID[24,TRIALHypoxia]"
[100] "r_FISHID[25,TRIALHypoxia]"
[101] "r_FISHID[26,TRIALHypoxia]"
[102] "r_FISHID[27,TRIALHypoxia]"
[103] "r_FISHID[28,TRIALHypoxia]"
[104] "r_FISHID[29,TRIALHypoxia]"
[105] "r_FISHID[30,TRIALHypoxia]"
[106] "r_FISHID[31,TRIALHypoxia]"
[107] "r_FISHID[32,TRIALHypoxia]"
[108] "r_FISHID[33,TRIALHypoxia]"
[109] "r_FISHID[34,TRIALHypoxia]"
[110] "r_FISHID[35,TRIALHypoxia]"
[111] "r_FISHID[36,TRIALHypoxia]"
[112] "r_FISHID[37,TRIALHypoxia]"
[113] "r_FISHID[38,TRIALHypoxia]"
[114] "r_FISHID[39,TRIALHypoxia]"
[115] "r_FISHID[40,TRIALHypoxia]"
[116] "r_FISHID[41,TRIALHypoxia]"
[117] "r_FISHID[42,TRIALHypoxia]"
[118] "r_FISHID[43,TRIALHypoxia]"
[119] "r_FISHID[44,TRIALHypoxia]"
[120] "r_FISHID[45,TRIALHypoxia]"
[121] "r_FISHID[46,TRIALHypoxia]"
[122] "r_FISHID[47,TRIALHypoxia]"
[123] "r_FISHID[48,TRIALHypoxia]"
[124] "r_FISHID[49,TRIALHypoxia]"
[125] "r_FISHID[50,TRIALHypoxia]"
[126] "r_FISHID[51,TRIALHypoxia]"
[127] "r_FISHID[52,TRIALHypoxia]"
[128] "r_FISHID[53,TRIALHypoxia]"
[129] "r_FISHID[54,TRIALHypoxia]"
[130] "r_FISHID[55,TRIALHypoxia]"
[131] "r_FISHID[56,TRIALHypoxia]"
[132] "r_FISHID[57,TRIALHypoxia]"
[133] "r_FISHID[58,TRIALHypoxia]"
[134] "r_FISHID[59,TRIALHypoxia]"
[135] "r_FISHID[60,TRIALHypoxia]"
[136] "r_FISHID[1,TRIALLowSalinity]"
[137] "r_FISHID[2,TRIALLowSalinity]"
[138] "r_FISHID[3,TRIALLowSalinity]"
[139] "r_FISHID[4,TRIALLowSalinity]"
[140] "r_FISHID[5,TRIALLowSalinity]"
[141] "r_FISHID[6,TRIALLowSalinity]"
[142] "r_FISHID[7,TRIALLowSalinity]"
[143] "r_FISHID[8,TRIALLowSalinity]"
[144] "r_FISHID[9,TRIALLowSalinity]"
[145] "r_FISHID[10,TRIALLowSalinity]"
[146] "r_FISHID[11,TRIALLowSalinity]"
[147] "r_FISHID[12,TRIALLowSalinity]"
[148] "r_FISHID[13,TRIALLowSalinity]"
[149] "r_FISHID[14,TRIALLowSalinity]"
[150] "r_FISHID[15,TRIALLowSalinity]"
[151] "r_FISHID[16,TRIALLowSalinity]"
[152] "r_FISHID[17,TRIALLowSalinity]"
[153] "r_FISHID[18,TRIALLowSalinity]"
[154] "r_FISHID[19,TRIALLowSalinity]"
[155] "r_FISHID[20,TRIALLowSalinity]"
[156] "r_FISHID[21,TRIALLowSalinity]"
[157] "r_FISHID[22,TRIALLowSalinity]"
[158] "r_FISHID[23,TRIALLowSalinity]"
[159] "r_FISHID[24,TRIALLowSalinity]"
[160] "r_FISHID[25,TRIALLowSalinity]"
[161] "r_FISHID[26,TRIALLowSalinity]"
[162] "r_FISHID[27,TRIALLowSalinity]"
[163] "r_FISHID[28,TRIALLowSalinity]"
[164] "r_FISHID[29,TRIALLowSalinity]"
[165] "r_FISHID[30,TRIALLowSalinity]"
[166] "r_FISHID[31,TRIALLowSalinity]"
[167] "r_FISHID[32,TRIALLowSalinity]"
[168] "r_FISHID[33,TRIALLowSalinity]"
[169] "r_FISHID[34,TRIALLowSalinity]"
[170] "r_FISHID[35,TRIALLowSalinity]"
[171] "r_FISHID[36,TRIALLowSalinity]"
[172] "r_FISHID[37,TRIALLowSalinity]"
[173] "r_FISHID[38,TRIALLowSalinity]"
[174] "r_FISHID[39,TRIALLowSalinity]"
[175] "r_FISHID[40,TRIALLowSalinity]"
[176] "r_FISHID[41,TRIALLowSalinity]"
[177] "r_FISHID[42,TRIALLowSalinity]"
[178] "r_FISHID[43,TRIALLowSalinity]"
[179] "r_FISHID[44,TRIALLowSalinity]"
[180] "r_FISHID[45,TRIALLowSalinity]"
[181] "r_FISHID[46,TRIALLowSalinity]"
[182] "r_FISHID[47,TRIALLowSalinity]"
[183] "r_FISHID[48,TRIALLowSalinity]"
[184] "r_FISHID[49,TRIALLowSalinity]"
[185] "r_FISHID[50,TRIALLowSalinity]"
[186] "r_FISHID[51,TRIALLowSalinity]"
[187] "r_FISHID[52,TRIALLowSalinity]"
[188] "r_FISHID[53,TRIALLowSalinity]"
[189] "r_FISHID[54,TRIALLowSalinity]"
[190] "r_FISHID[55,TRIALLowSalinity]"
[191] "r_FISHID[56,TRIALLowSalinity]"
[192] "r_FISHID[57,TRIALLowSalinity]"
[193] "r_FISHID[58,TRIALLowSalinity]"
[194] "r_FISHID[59,TRIALLowSalinity]"
[195] "r_FISHID[60,TRIALLowSalinity]"
[196] "prior_Intercept"
[197] "prior_b_TRIALHypoxia"
[198] "prior_b_TRIALLowSalinity"
[199] "prior_b_SMR_contr"
[200] "prior_b_MASS"
[201] "prior_sigma"
[202] "prior_sd_FISHID"
[203] "prior_cor_FISHID"
[204] "lprior"
[205] "lp__"
[206] "accept_stat__"
[207] "stepsize__"
[208] "treedepth__"
[209] "n_leapfrog__"
[210] "divergent__"
[211] "energy__"
norin.brm4 |>
posterior_samples() |>
dplyr::select(-`lp__`) |>
pivot_longer(everything(), names_to = "key") |>
filter(!str_detect(key, "^r")) |>
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
Class = case_when(
str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
str_detect(key, "b_TRIAL.*|prior_b_TRIAL.*") & !str_detect(key, ".*\\:.*") ~ "TRIAL",
str_detect(key, "b_SMR_contr|prior_b_SMR_contr") ~ "SMR",
str_detect(key, "b_MASS|prior_b_MASS") ~ "MASS",
str_detect(key, ".*\\:.*|prior_b_.*\\:.*") ~ "Interaction",
str_detect(key, "sd") ~ "sd",
str_detect(key, "^cor|prior_cor") ~ "cor",
str_detect(key, "sigma") ~ "sigma"
),
Par = str_replace(key, "b_", "")
) |>
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge()) +
facet_wrap(~Class, scales = "free")Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
We can compare the two models using LOO
Computed from 2400 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -832.1 13.7
p_loo 24.5 3.9
looic 1664.2 27.4
------
MCSE of elpd_loo is 0.2.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.
Warning: Found 54 observations with a pareto_k > 0.7 in model 'norin.brm4'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.
Computed from 2400 by 180 log-likelihood matrix.
Estimate SE
elpd_loo -780.1 9.2
p_loo 85.7 6.1
looic 1560.3 18.4
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 0.7]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 126 70.0% 47
(0.7, 1] (bad) 47 26.1% <NA>
(1, Inf) (very bad) 7 3.9% <NA>
See help('pareto-k-diagnostic') for details.
elpd_diff se_diff
norin.brm4 0.0 0.0
norin.brm3 -52.0 10.2
There is substantially more support for the more complex random intercept/slope model over the simpler random intercept only model. Consequently, we will continue with the random intercept/slope model.
6 MCMC sampling diagnostics
The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.
See list of available diagnostics by name
bayesplot MCMC module:
mcmc_acf
mcmc_acf_bar
mcmc_areas
mcmc_areas_ridges
mcmc_combo
mcmc_dens
mcmc_dens_chains
mcmc_dens_overlay
mcmc_hex
mcmc_hist
mcmc_hist_by_chain
mcmc_intervals
mcmc_neff
mcmc_neff_hist
mcmc_nuts_acceptance
mcmc_nuts_divergence
mcmc_nuts_energy
mcmc_nuts_stepsize
mcmc_nuts_treedepth
mcmc_pairs
mcmc_parcoord
mcmc_rank_ecdf
mcmc_rank_hist
mcmc_rank_overlay
mcmc_recover_hist
mcmc_recover_intervals
mcmc_recover_scatter
mcmc_rhat
mcmc_rhat_hist
mcmc_scatter
mcmc_trace
mcmc_trace_highlight
mcmc_violin
Of these, we will focus on:
- trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
The chains appear well mixed and very similar
- acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
- rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
All Rhat values are below 1.05, suggesting the chains have converged.
neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
[1] "b_Intercept"
[2] "b_TRIALHypoxia"
[3] "b_TRIALLowSalinity"
[4] "b_SMR_contr"
[5] "b_MASS"
[6] "b_TRIALHypoxia:SMR_contr"
[7] "b_TRIALLowSalinity:SMR_contr"
[8] "sd_FISHID__Intercept"
[9] "sd_FISHID__TRIALHypoxia"
[10] "sd_FISHID__TRIALLowSalinity"
[11] "cor_FISHID__Intercept__TRIALHypoxia"
[12] "cor_FISHID__Intercept__TRIALLowSalinity"
[13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
[14] "sigma"
[15] "Intercept"
[16] "r_FISHID[1,Intercept]"
[17] "r_FISHID[2,Intercept]"
[18] "r_FISHID[3,Intercept]"
[19] "r_FISHID[4,Intercept]"
[20] "r_FISHID[5,Intercept]"
[21] "r_FISHID[6,Intercept]"
[22] "r_FISHID[7,Intercept]"
[23] "r_FISHID[8,Intercept]"
[24] "r_FISHID[9,Intercept]"
[25] "r_FISHID[10,Intercept]"
[26] "r_FISHID[11,Intercept]"
[27] "r_FISHID[12,Intercept]"
[28] "r_FISHID[13,Intercept]"
[29] "r_FISHID[14,Intercept]"
[30] "r_FISHID[15,Intercept]"
[31] "r_FISHID[16,Intercept]"
[32] "r_FISHID[17,Intercept]"
[33] "r_FISHID[18,Intercept]"
[34] "r_FISHID[19,Intercept]"
[35] "r_FISHID[20,Intercept]"
[36] "r_FISHID[21,Intercept]"
[37] "r_FISHID[22,Intercept]"
[38] "r_FISHID[23,Intercept]"
[39] "r_FISHID[24,Intercept]"
[40] "r_FISHID[25,Intercept]"
[41] "r_FISHID[26,Intercept]"
[42] "r_FISHID[27,Intercept]"
[43] "r_FISHID[28,Intercept]"
[44] "r_FISHID[29,Intercept]"
[45] "r_FISHID[30,Intercept]"
[46] "r_FISHID[31,Intercept]"
[47] "r_FISHID[32,Intercept]"
[48] "r_FISHID[33,Intercept]"
[49] "r_FISHID[34,Intercept]"
[50] "r_FISHID[35,Intercept]"
[51] "r_FISHID[36,Intercept]"
[52] "r_FISHID[37,Intercept]"
[53] "r_FISHID[38,Intercept]"
[54] "r_FISHID[39,Intercept]"
[55] "r_FISHID[40,Intercept]"
[56] "r_FISHID[41,Intercept]"
[57] "r_FISHID[42,Intercept]"
[58] "r_FISHID[43,Intercept]"
[59] "r_FISHID[44,Intercept]"
[60] "r_FISHID[45,Intercept]"
[61] "r_FISHID[46,Intercept]"
[62] "r_FISHID[47,Intercept]"
[63] "r_FISHID[48,Intercept]"
[64] "r_FISHID[49,Intercept]"
[65] "r_FISHID[50,Intercept]"
[66] "r_FISHID[51,Intercept]"
[67] "r_FISHID[52,Intercept]"
[68] "r_FISHID[53,Intercept]"
[69] "r_FISHID[54,Intercept]"
[70] "r_FISHID[55,Intercept]"
[71] "r_FISHID[56,Intercept]"
[72] "r_FISHID[57,Intercept]"
[73] "r_FISHID[58,Intercept]"
[74] "r_FISHID[59,Intercept]"
[75] "r_FISHID[60,Intercept]"
[76] "r_FISHID[1,TRIALHypoxia]"
[77] "r_FISHID[2,TRIALHypoxia]"
[78] "r_FISHID[3,TRIALHypoxia]"
[79] "r_FISHID[4,TRIALHypoxia]"
[80] "r_FISHID[5,TRIALHypoxia]"
[81] "r_FISHID[6,TRIALHypoxia]"
[82] "r_FISHID[7,TRIALHypoxia]"
[83] "r_FISHID[8,TRIALHypoxia]"
[84] "r_FISHID[9,TRIALHypoxia]"
[85] "r_FISHID[10,TRIALHypoxia]"
[86] "r_FISHID[11,TRIALHypoxia]"
[87] "r_FISHID[12,TRIALHypoxia]"
[88] "r_FISHID[13,TRIALHypoxia]"
[89] "r_FISHID[14,TRIALHypoxia]"
[90] "r_FISHID[15,TRIALHypoxia]"
[91] "r_FISHID[16,TRIALHypoxia]"
[92] "r_FISHID[17,TRIALHypoxia]"
[93] "r_FISHID[18,TRIALHypoxia]"
[94] "r_FISHID[19,TRIALHypoxia]"
[95] "r_FISHID[20,TRIALHypoxia]"
[96] "r_FISHID[21,TRIALHypoxia]"
[97] "r_FISHID[22,TRIALHypoxia]"
[98] "r_FISHID[23,TRIALHypoxia]"
[99] "r_FISHID[24,TRIALHypoxia]"
[100] "r_FISHID[25,TRIALHypoxia]"
[101] "r_FISHID[26,TRIALHypoxia]"
[102] "r_FISHID[27,TRIALHypoxia]"
[103] "r_FISHID[28,TRIALHypoxia]"
[104] "r_FISHID[29,TRIALHypoxia]"
[105] "r_FISHID[30,TRIALHypoxia]"
[106] "r_FISHID[31,TRIALHypoxia]"
[107] "r_FISHID[32,TRIALHypoxia]"
[108] "r_FISHID[33,TRIALHypoxia]"
[109] "r_FISHID[34,TRIALHypoxia]"
[110] "r_FISHID[35,TRIALHypoxia]"
[111] "r_FISHID[36,TRIALHypoxia]"
[112] "r_FISHID[37,TRIALHypoxia]"
[113] "r_FISHID[38,TRIALHypoxia]"
[114] "r_FISHID[39,TRIALHypoxia]"
[115] "r_FISHID[40,TRIALHypoxia]"
[116] "r_FISHID[41,TRIALHypoxia]"
[117] "r_FISHID[42,TRIALHypoxia]"
[118] "r_FISHID[43,TRIALHypoxia]"
[119] "r_FISHID[44,TRIALHypoxia]"
[120] "r_FISHID[45,TRIALHypoxia]"
[121] "r_FISHID[46,TRIALHypoxia]"
[122] "r_FISHID[47,TRIALHypoxia]"
[123] "r_FISHID[48,TRIALHypoxia]"
[124] "r_FISHID[49,TRIALHypoxia]"
[125] "r_FISHID[50,TRIALHypoxia]"
[126] "r_FISHID[51,TRIALHypoxia]"
[127] "r_FISHID[52,TRIALHypoxia]"
[128] "r_FISHID[53,TRIALHypoxia]"
[129] "r_FISHID[54,TRIALHypoxia]"
[130] "r_FISHID[55,TRIALHypoxia]"
[131] "r_FISHID[56,TRIALHypoxia]"
[132] "r_FISHID[57,TRIALHypoxia]"
[133] "r_FISHID[58,TRIALHypoxia]"
[134] "r_FISHID[59,TRIALHypoxia]"
[135] "r_FISHID[60,TRIALHypoxia]"
[136] "r_FISHID[1,TRIALLowSalinity]"
[137] "r_FISHID[2,TRIALLowSalinity]"
[138] "r_FISHID[3,TRIALLowSalinity]"
[139] "r_FISHID[4,TRIALLowSalinity]"
[140] "r_FISHID[5,TRIALLowSalinity]"
[141] "r_FISHID[6,TRIALLowSalinity]"
[142] "r_FISHID[7,TRIALLowSalinity]"
[143] "r_FISHID[8,TRIALLowSalinity]"
[144] "r_FISHID[9,TRIALLowSalinity]"
[145] "r_FISHID[10,TRIALLowSalinity]"
[146] "r_FISHID[11,TRIALLowSalinity]"
[147] "r_FISHID[12,TRIALLowSalinity]"
[148] "r_FISHID[13,TRIALLowSalinity]"
[149] "r_FISHID[14,TRIALLowSalinity]"
[150] "r_FISHID[15,TRIALLowSalinity]"
[151] "r_FISHID[16,TRIALLowSalinity]"
[152] "r_FISHID[17,TRIALLowSalinity]"
[153] "r_FISHID[18,TRIALLowSalinity]"
[154] "r_FISHID[19,TRIALLowSalinity]"
[155] "r_FISHID[20,TRIALLowSalinity]"
[156] "r_FISHID[21,TRIALLowSalinity]"
[157] "r_FISHID[22,TRIALLowSalinity]"
[158] "r_FISHID[23,TRIALLowSalinity]"
[159] "r_FISHID[24,TRIALLowSalinity]"
[160] "r_FISHID[25,TRIALLowSalinity]"
[161] "r_FISHID[26,TRIALLowSalinity]"
[162] "r_FISHID[27,TRIALLowSalinity]"
[163] "r_FISHID[28,TRIALLowSalinity]"
[164] "r_FISHID[29,TRIALLowSalinity]"
[165] "r_FISHID[30,TRIALLowSalinity]"
[166] "r_FISHID[31,TRIALLowSalinity]"
[167] "r_FISHID[32,TRIALLowSalinity]"
[168] "r_FISHID[33,TRIALLowSalinity]"
[169] "r_FISHID[34,TRIALLowSalinity]"
[170] "r_FISHID[35,TRIALLowSalinity]"
[171] "r_FISHID[36,TRIALLowSalinity]"
[172] "r_FISHID[37,TRIALLowSalinity]"
[173] "r_FISHID[38,TRIALLowSalinity]"
[174] "r_FISHID[39,TRIALLowSalinity]"
[175] "r_FISHID[40,TRIALLowSalinity]"
[176] "r_FISHID[41,TRIALLowSalinity]"
[177] "r_FISHID[42,TRIALLowSalinity]"
[178] "r_FISHID[43,TRIALLowSalinity]"
[179] "r_FISHID[44,TRIALLowSalinity]"
[180] "r_FISHID[45,TRIALLowSalinity]"
[181] "r_FISHID[46,TRIALLowSalinity]"
[182] "r_FISHID[47,TRIALLowSalinity]"
[183] "r_FISHID[48,TRIALLowSalinity]"
[184] "r_FISHID[49,TRIALLowSalinity]"
[185] "r_FISHID[50,TRIALLowSalinity]"
[186] "r_FISHID[51,TRIALLowSalinity]"
[187] "r_FISHID[52,TRIALLowSalinity]"
[188] "r_FISHID[53,TRIALLowSalinity]"
[189] "r_FISHID[54,TRIALLowSalinity]"
[190] "r_FISHID[55,TRIALLowSalinity]"
[191] "r_FISHID[56,TRIALLowSalinity]"
[192] "r_FISHID[57,TRIALLowSalinity]"
[193] "r_FISHID[58,TRIALLowSalinity]"
[194] "r_FISHID[59,TRIALLowSalinity]"
[195] "r_FISHID[60,TRIALLowSalinity]"
[196] "prior_Intercept"
[197] "prior_b_TRIALHypoxia"
[198] "prior_b_TRIALLowSalinity"
[199] "prior_b_SMR_contr"
[200] "prior_b_MASS"
[201] "prior_sigma"
[202] "prior_sd_FISHID"
[203] "prior_cor_FISHID"
[204] "lprior"
[205] "lp__"
[206] "accept_stat__"
[207] "stepsize__"
[208] "treedepth__"
[209] "n_leapfrog__"
[210] "divergent__"
[211] "energy__"
pars <- norin.brm4 |> get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") |> na.omit()
norin.brm4$fit |>
stan_trace(pars = pars)The chains appear well mixed and very similar
- stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.05, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.
The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).
Warning in custom.sort(D$Parameter): NAs introduced by coercion
The chains appear well mixed and very similar
- gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
There is no evidence of autocorrelation in the MCMC samples
- stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.
All Rhat values are below 1.01, suggesting the chains have converged.
stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).
If the ratios are low, tightening the priors may help.
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
dplyr 1.1.0.
ℹ Please use `reframe()` instead.
ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
always returns an ungrouped data frame and adjust accordingly.
ℹ The deprecated feature was likely used in the ggmcmc package.
Please report the issue at <https://github.com/xfim/ggmcmc/issues/>.
Ratios all very high.
Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.
See list of available diagnostics by name
bayesplot PPC module:
ppc_bars
ppc_bars_grouped
ppc_boxplot
ppc_dens
ppc_dens_overlay
ppc_dens_overlay_grouped
ppc_ecdf_overlay
ppc_ecdf_overlay_grouped
ppc_error_binned
ppc_error_hist
ppc_error_hist_grouped
ppc_error_scatter
ppc_error_scatter_avg
ppc_error_scatter_avg_grouped
ppc_error_scatter_avg_vs_x
ppc_freqpoly
ppc_freqpoly_grouped
ppc_hist
ppc_intervals
ppc_intervals_grouped
ppc_km_overlay
ppc_km_overlay_grouped
ppc_loo_intervals
ppc_loo_pit_ecdf
ppc_loo_pit_overlay
ppc_loo_pit_qq
ppc_loo_ribbon
ppc_pit_ecdf
ppc_pit_ecdf_grouped
ppc_ribbon
ppc_ribbon_grouped
ppc_rootogram
ppc_scatter
ppc_scatter_avg
ppc_scatter_avg_grouped
ppc_stat
ppc_stat_2d
ppc_stat_freqpoly
ppc_stat_freqpoly_grouped
ppc_stat_grouped
ppc_violin_grouped
- dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.
The model draws appear to be consistent with the observed data.
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response
Using all posterior draws for ppc type 'error_scatter_avg' by default.
This is not really interpretable
- intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
Using all posterior draws for ppc type 'intervals' by default.
Warning: The following arguments were unrecognized and ignored: group
Using all posterior draws for ppc type 'intervals_grouped' by default.
Using all posterior draws for ppc type 'violin_grouped' by default.
The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.
DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.
We need to supply:
- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation
preds <- norin.brm4 |> posterior_predict(ndraws = 250, summary = FALSE)
norin.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE
)
plot(norin.resids, quantreg = FALSE)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.2998, p-value < 2.2e-16
alternative hypothesis: two.sided
norin.resids <- make_brms_dharma_res(norin.brm4, integerResponse = FALSE)
wrap_elements(~ testUniformity(norin.resids)) +
wrap_elements(~ plotResiduals(norin.resids, form = factor(rep(1, nrow(norin))))) +
wrap_elements(~ plotResiduals(norin.resids)) +
wrap_elements(~ testDispersion(norin.resids))Conclusions:
- the simulated residuals do suggest some issues with the residuals
- the model appears to be under-dispersed
- there is evidence of a lack of fit.
7 Partial effects plots
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
It is not really possible to do this via the fitted draws as it would not be marginalising over MASS or FISHID.
8 Model investigation
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
Family: gaussian
Links: mu = identity; sigma = identity
Formula: CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS
Data: norin (Number of observations: 180)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~FISHID (Number of levels: 60)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept) 11.77 2.94 5.72 17.45 1.01
sd(TRIALHypoxia) 6.31 3.72 0.39 14.59 1.00
sd(TRIALLowSalinity) 28.09 3.90 20.80 36.43 1.00
cor(Intercept,TRIALHypoxia) -0.33 0.38 -0.88 0.62 1.00
cor(Intercept,TRIALLowSalinity) 0.08 0.24 -0.34 0.61 1.00
cor(TRIALHypoxia,TRIALLowSalinity) -0.48 0.37 -0.95 0.52 1.00
Bulk_ESS Tail_ESS
sd(Intercept) 729 1248
sd(TRIALHypoxia) 725 672
sd(TRIALLowSalinity) 1501 1742
cor(Intercept,TRIALHypoxia) 1857 2050
cor(Intercept,TRIALLowSalinity) 1368 1600
cor(TRIALHypoxia,TRIALLowSalinity) 611 519
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 147.10 22.28 103.34 190.74 1.00 2420
TRIALHypoxia -107.54 21.45 -149.71 -64.56 1.00 2372
TRIALLowSalinity -72.63 32.21 -134.70 -5.92 1.00 2262
SMR_contr -19.59 3.68 -26.60 -12.29 1.00 2433
MASS 0.12 0.24 -0.34 0.59 1.00 2409
TRIALHypoxia:SMR_contr 10.96 4.14 2.86 18.97 1.00 2370
TRIALLowSalinity:SMR_contr 5.95 6.19 -6.63 17.75 1.00 2278
Tail_ESS
Intercept 2292
TRIALHypoxia 2227
TRIALLowSalinity 2002
SMR_contr 1999
MASS 2288
TRIALHypoxia:SMR_contr 2383
TRIALLowSalinity:SMR_contr 2159
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 13.37 1.94 9.50 17.21 1.01 549 517
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Conclusions:
the intercept indicates that the average metabolic change associated with the High Temperature trial (at the mean SMR_contr and MASS) across all fish was 147.1.
for the same SMR_contr and MASS, the average metabolic change was 107.54 less in the Hypoxia trial and 72.63 less in the Salinity trial than the High Temperature trial.
in the High Temperature trial, every 1 unit increase in SMR_contr was associated with a decline in metabolic change of
round(-1*norin.sum$fixed[4,1], 2)units.in the High Temperature trial (and mean SMR_contr), every 1 unit increase in MASS was associated with an increase in metabolic change of
round(-1*norin.sum$fixed[5,1], 2)units.the rate of change in metabolic change associated with the Hypoxia trial was
round(-1*norin.sum$fixed[6,1], 2)units higher than that of the High Temperature trial andround(-1*norin.sum$fixed[7,1], 2)units higher in the Low Salinity trial.the variance in intercepts across all fish is 11.77
the variance in responses between High Temperature and Hypoxia trials across all fish is 6.31.
the variance in responses between High Temperature and Low Salinity trials across all fish is 28.09.
the scale of variance between fish within a treatment
(sigma, 13.37) is substantially higher than that both within the High Temperature trial and between the High Temperature and Hypoxia trials, yet lower than that between High Temperature and Low Salinity trials.
# A tibble: 205 × 7
variable median lower upper rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 146. 1.04e+2 191. 1.00 2420. 2292.
2 b_TRIALHypoxia -108. -1.48e+2 -64.1 1.00 2372. 2227.
3 b_TRIALLowSalinity -72.9 -1.37e+2 -10.2 1.000 2262. 2002.
4 b_SMR_contr -19.7 -2.67e+1 -12.4 1.000 2433. 1999.
5 b_MASS 0.118 -3.51e-1 0.573 1.000 2408. 2288.
6 b_TRIALHypoxia:SMR_contr 11.0 2.84e+0 18.9 1.00 2370. 2383.
7 b_TRIALLowSalinity:SMR_con… 5.89 -6.37e+0 17.8 0.999 2279. 2159.
8 sd_FISHID__Intercept 11.8 5.64e+0 17.3 1.01 729. 1248.
9 sd_FISHID__TRIALHypoxia 6.07 2.69e-3 12.8 1.00 725. 672.
10 sd_FISHID__TRIALLowSalinity 28.0 2.03e+1 35.5 1.00 1500. 1742.
# ℹ 195 more rows
## or if you want to exclude some parameters
norin.brm4 |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk,
ess_tail
) |>
filter(str_detect(variable, "prior|^r_|^lp__", negate = TRUE))# A tibble: 15 × 7
variable median lower upper rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 1.46e+2 1.04e+2 191. 1.00 2420. 2292.
2 b_TRIALHypoxia -1.08e+2 -1.48e+2 -64.1 1.00 2372. 2227.
3 b_TRIALLowSalinity -7.29e+1 -1.37e+2 -10.2 1.000 2262. 2002.
4 b_SMR_contr -1.97e+1 -2.67e+1 -12.4 1.000 2433. 1999.
5 b_MASS 1.18e-1 -3.51e-1 0.573 1.000 2408. 2288.
6 b_TRIALHypoxia:SMR_contr 1.10e+1 2.84e+0 18.9 1.00 2370. 2383.
7 b_TRIALLowSalinity:SMR_con… 5.89e+0 -6.37e+0 17.8 0.999 2279. 2159.
8 sd_FISHID__Intercept 1.18e+1 5.64e+0 17.3 1.01 729. 1248.
9 sd_FISHID__TRIALHypoxia 6.07e+0 2.69e-3 12.8 1.00 725. 672.
10 sd_FISHID__TRIALLowSalinity 2.80e+1 2.03e+1 35.5 1.00 1500. 1742.
11 cor_FISHID__Intercept__TRI… -3.97e-1 -9.58e-1 0.446 1.00 1856. 2050.
12 cor_FISHID__Intercept__TRI… 4.94e-2 -3.61e-1 0.574 1.00 1369. 1600.
13 cor_FISHID__TRIALHypoxia__… -5.50e-1 -9.99e-1 0.249 1.00 611. 519.
14 sigma 1.34e+1 9.83e+0 17.4 1.01 548. 517.
15 Intercept 1.97e+1 1.54e+1 23.5 0.999 2504. 2228.
# A draws_df: 800 iterations, 3 chains, and 205 variables
b_Intercept b_TRIALHypoxia b_TRIALLowSalinity b_SMR_contr b_MASS
1 139 -153 -73 -20 0.3392
2 132 -103 -73 -17 0.2897
3 135 -107 -39 -19 0.3357
4 169 -140 -81 -23 -0.0097
5 144 -127 -3 -19 0.1346
6 154 -104 1 -22 0.2451
7 178 -105 -48 -24 -0.0865
8 155 -105 -56 -21 0.0615
9 132 -80 -62 -15 0.0266
10 149 -106 -109 -21 0.2686
b_TRIALHypoxia:SMR_contr b_TRIALLowSalinity:SMR_contr sd_FISHID__Intercept
1 19.9 7.1 10.8
2 9.8 6.5 10.8
3 10.0 -1.2 10.0
4 17.0 7.6 15.0
5 14.9 -7.3 12.0
6 10.3 -9.0 11.5
7 10.3 1.8 10.3
8 10.7 3.1 13.0
9 4.8 4.2 15.1
10 10.3 12.5 8.3
# ... with 2390 more draws, and 197 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
norin.brm4 |>
as_draws_df() |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk,
ess_tail
)# A tibble: 205 × 7
variable median lower upper rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 146. 1.04e+2 191. 1.00 2420. 2292.
2 b_TRIALHypoxia -108. -1.48e+2 -64.1 1.00 2372. 2227.
3 b_TRIALLowSalinity -72.9 -1.37e+2 -10.2 1.000 2262. 2002.
4 b_SMR_contr -19.7 -2.67e+1 -12.4 1.000 2433. 1999.
5 b_MASS 0.118 -3.51e-1 0.573 1.000 2408. 2288.
6 b_TRIALHypoxia:SMR_contr 11.0 2.84e+0 18.9 1.00 2370. 2383.
7 b_TRIALLowSalinity:SMR_con… 5.89 -6.37e+0 17.8 0.999 2279. 2159.
8 sd_FISHID__Intercept 11.8 5.64e+0 17.3 1.01 729. 1248.
9 sd_FISHID__TRIALHypoxia 6.07 2.69e-3 12.8 1.00 725. 672.
10 sd_FISHID__TRIALLowSalinity 28.0 2.03e+1 35.5 1.00 1500. 1742.
# ℹ 195 more rows
## or if you want to exclude some parameters
norin.brm4 |>
as_draws_df() |>
summarise_draws(
median,
~ HDInterval::hdi(.x),
rhat,
ess_bulk,
ess_tail
) |>
filter(str_detect(variable, "prior|^r_|^lp__", negate = TRUE))# A tibble: 15 × 7
variable median lower upper rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 b_Intercept 1.46e+2 1.04e+2 191. 1.00 2420. 2292.
2 b_TRIALHypoxia -1.08e+2 -1.48e+2 -64.1 1.00 2372. 2227.
3 b_TRIALLowSalinity -7.29e+1 -1.37e+2 -10.2 1.000 2262. 2002.
4 b_SMR_contr -1.97e+1 -2.67e+1 -12.4 1.000 2433. 1999.
5 b_MASS 1.18e-1 -3.51e-1 0.573 1.000 2408. 2288.
6 b_TRIALHypoxia:SMR_contr 1.10e+1 2.84e+0 18.9 1.00 2370. 2383.
7 b_TRIALLowSalinity:SMR_con… 5.89e+0 -6.37e+0 17.8 0.999 2279. 2159.
8 sd_FISHID__Intercept 1.18e+1 5.64e+0 17.3 1.01 729. 1248.
9 sd_FISHID__TRIALHypoxia 6.07e+0 2.69e-3 12.8 1.00 725. 672.
10 sd_FISHID__TRIALLowSalinity 2.80e+1 2.03e+1 35.5 1.00 1500. 1742.
11 cor_FISHID__Intercept__TRI… -3.97e-1 -9.58e-1 0.446 1.00 1856. 2050.
12 cor_FISHID__Intercept__TRI… 4.94e-2 -3.61e-1 0.574 1.00 1369. 1600.
13 cor_FISHID__TRIALHypoxia__… -5.50e-1 -9.99e-1 0.249 1.00 611. 519.
14 sigma 1.34e+1 9.83e+0 17.4 1.01 548. 517.
15 Intercept 1.97e+1 1.54e+1 23.5 0.999 2504. 2228.
norin.brm4$fit |>
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE
)# A tibble: 204 × 7
term estimate std.error conf.low conf.high rhat ess
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 b_Intercept 147. 22.3 1.04e+2 191. 0.999 2412
2 b_TRIALHypoxia -108. 21.4 -1.48e+2 -64.1 1.000 2366
3 b_TRIALLowSalinity -72.6 32.2 -1.37e+2 -10.2 0.999 2234
4 b_SMR_contr -19.6 3.68 -2.67e+1 -12.4 0.999 2427
5 b_MASS 0.121 0.240 -3.51e-1 0.573 1.000 2414
6 b_TRIALHypoxia:SMR_contr 11.0 4.14 2.84e+0 18.9 1.000 2358
7 b_TRIALLowSalinity:SMR_con… 5.95 6.19 -6.37e+0 17.8 0.999 2245
8 sd_FISHID__Intercept 11.8 2.94 5.64e+0 17.3 1.01 633
9 sd_FISHID__TRIALHypoxia 6.31 3.72 2.69e-3 12.8 1.00 704
10 sd_FISHID__TRIALLowSalinity 28.1 3.90 2.03e+1 35.5 1.00 1499
# ℹ 194 more rows
Conclusions:
see above
[1] "b_Intercept"
[2] "b_TRIALHypoxia"
[3] "b_TRIALLowSalinity"
[4] "b_SMR_contr"
[5] "b_MASS"
[6] "b_TRIALHypoxia:SMR_contr"
[7] "b_TRIALLowSalinity:SMR_contr"
[8] "sd_FISHID__Intercept"
[9] "sd_FISHID__TRIALHypoxia"
[10] "sd_FISHID__TRIALLowSalinity"
[11] "cor_FISHID__Intercept__TRIALHypoxia"
[12] "cor_FISHID__Intercept__TRIALLowSalinity"
[13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
[14] "sigma"
[15] "Intercept"
[16] "r_FISHID[1,Intercept]"
[17] "r_FISHID[2,Intercept]"
[18] "r_FISHID[3,Intercept]"
[19] "r_FISHID[4,Intercept]"
[20] "r_FISHID[5,Intercept]"
[21] "r_FISHID[6,Intercept]"
[22] "r_FISHID[7,Intercept]"
[23] "r_FISHID[8,Intercept]"
[24] "r_FISHID[9,Intercept]"
[25] "r_FISHID[10,Intercept]"
[26] "r_FISHID[11,Intercept]"
[27] "r_FISHID[12,Intercept]"
[28] "r_FISHID[13,Intercept]"
[29] "r_FISHID[14,Intercept]"
[30] "r_FISHID[15,Intercept]"
[31] "r_FISHID[16,Intercept]"
[32] "r_FISHID[17,Intercept]"
[33] "r_FISHID[18,Intercept]"
[34] "r_FISHID[19,Intercept]"
[35] "r_FISHID[20,Intercept]"
[36] "r_FISHID[21,Intercept]"
[37] "r_FISHID[22,Intercept]"
[38] "r_FISHID[23,Intercept]"
[39] "r_FISHID[24,Intercept]"
[40] "r_FISHID[25,Intercept]"
[41] "r_FISHID[26,Intercept]"
[42] "r_FISHID[27,Intercept]"
[43] "r_FISHID[28,Intercept]"
[44] "r_FISHID[29,Intercept]"
[45] "r_FISHID[30,Intercept]"
[46] "r_FISHID[31,Intercept]"
[47] "r_FISHID[32,Intercept]"
[48] "r_FISHID[33,Intercept]"
[49] "r_FISHID[34,Intercept]"
[50] "r_FISHID[35,Intercept]"
[51] "r_FISHID[36,Intercept]"
[52] "r_FISHID[37,Intercept]"
[53] "r_FISHID[38,Intercept]"
[54] "r_FISHID[39,Intercept]"
[55] "r_FISHID[40,Intercept]"
[56] "r_FISHID[41,Intercept]"
[57] "r_FISHID[42,Intercept]"
[58] "r_FISHID[43,Intercept]"
[59] "r_FISHID[44,Intercept]"
[60] "r_FISHID[45,Intercept]"
[61] "r_FISHID[46,Intercept]"
[62] "r_FISHID[47,Intercept]"
[63] "r_FISHID[48,Intercept]"
[64] "r_FISHID[49,Intercept]"
[65] "r_FISHID[50,Intercept]"
[66] "r_FISHID[51,Intercept]"
[67] "r_FISHID[52,Intercept]"
[68] "r_FISHID[53,Intercept]"
[69] "r_FISHID[54,Intercept]"
[70] "r_FISHID[55,Intercept]"
[71] "r_FISHID[56,Intercept]"
[72] "r_FISHID[57,Intercept]"
[73] "r_FISHID[58,Intercept]"
[74] "r_FISHID[59,Intercept]"
[75] "r_FISHID[60,Intercept]"
[76] "r_FISHID[1,TRIALHypoxia]"
[77] "r_FISHID[2,TRIALHypoxia]"
[78] "r_FISHID[3,TRIALHypoxia]"
[79] "r_FISHID[4,TRIALHypoxia]"
[80] "r_FISHID[5,TRIALHypoxia]"
[81] "r_FISHID[6,TRIALHypoxia]"
[82] "r_FISHID[7,TRIALHypoxia]"
[83] "r_FISHID[8,TRIALHypoxia]"
[84] "r_FISHID[9,TRIALHypoxia]"
[85] "r_FISHID[10,TRIALHypoxia]"
[86] "r_FISHID[11,TRIALHypoxia]"
[87] "r_FISHID[12,TRIALHypoxia]"
[88] "r_FISHID[13,TRIALHypoxia]"
[89] "r_FISHID[14,TRIALHypoxia]"
[90] "r_FISHID[15,TRIALHypoxia]"
[91] "r_FISHID[16,TRIALHypoxia]"
[92] "r_FISHID[17,TRIALHypoxia]"
[93] "r_FISHID[18,TRIALHypoxia]"
[94] "r_FISHID[19,TRIALHypoxia]"
[95] "r_FISHID[20,TRIALHypoxia]"
[96] "r_FISHID[21,TRIALHypoxia]"
[97] "r_FISHID[22,TRIALHypoxia]"
[98] "r_FISHID[23,TRIALHypoxia]"
[99] "r_FISHID[24,TRIALHypoxia]"
[100] "r_FISHID[25,TRIALHypoxia]"
[101] "r_FISHID[26,TRIALHypoxia]"
[102] "r_FISHID[27,TRIALHypoxia]"
[103] "r_FISHID[28,TRIALHypoxia]"
[104] "r_FISHID[29,TRIALHypoxia]"
[105] "r_FISHID[30,TRIALHypoxia]"
[106] "r_FISHID[31,TRIALHypoxia]"
[107] "r_FISHID[32,TRIALHypoxia]"
[108] "r_FISHID[33,TRIALHypoxia]"
[109] "r_FISHID[34,TRIALHypoxia]"
[110] "r_FISHID[35,TRIALHypoxia]"
[111] "r_FISHID[36,TRIALHypoxia]"
[112] "r_FISHID[37,TRIALHypoxia]"
[113] "r_FISHID[38,TRIALHypoxia]"
[114] "r_FISHID[39,TRIALHypoxia]"
[115] "r_FISHID[40,TRIALHypoxia]"
[116] "r_FISHID[41,TRIALHypoxia]"
[117] "r_FISHID[42,TRIALHypoxia]"
[118] "r_FISHID[43,TRIALHypoxia]"
[119] "r_FISHID[44,TRIALHypoxia]"
[120] "r_FISHID[45,TRIALHypoxia]"
[121] "r_FISHID[46,TRIALHypoxia]"
[122] "r_FISHID[47,TRIALHypoxia]"
[123] "r_FISHID[48,TRIALHypoxia]"
[124] "r_FISHID[49,TRIALHypoxia]"
[125] "r_FISHID[50,TRIALHypoxia]"
[126] "r_FISHID[51,TRIALHypoxia]"
[127] "r_FISHID[52,TRIALHypoxia]"
[128] "r_FISHID[53,TRIALHypoxia]"
[129] "r_FISHID[54,TRIALHypoxia]"
[130] "r_FISHID[55,TRIALHypoxia]"
[131] "r_FISHID[56,TRIALHypoxia]"
[132] "r_FISHID[57,TRIALHypoxia]"
[133] "r_FISHID[58,TRIALHypoxia]"
[134] "r_FISHID[59,TRIALHypoxia]"
[135] "r_FISHID[60,TRIALHypoxia]"
[136] "r_FISHID[1,TRIALLowSalinity]"
[137] "r_FISHID[2,TRIALLowSalinity]"
[138] "r_FISHID[3,TRIALLowSalinity]"
[139] "r_FISHID[4,TRIALLowSalinity]"
[140] "r_FISHID[5,TRIALLowSalinity]"
[141] "r_FISHID[6,TRIALLowSalinity]"
[142] "r_FISHID[7,TRIALLowSalinity]"
[143] "r_FISHID[8,TRIALLowSalinity]"
[144] "r_FISHID[9,TRIALLowSalinity]"
[145] "r_FISHID[10,TRIALLowSalinity]"
[146] "r_FISHID[11,TRIALLowSalinity]"
[147] "r_FISHID[12,TRIALLowSalinity]"
[148] "r_FISHID[13,TRIALLowSalinity]"
[149] "r_FISHID[14,TRIALLowSalinity]"
[150] "r_FISHID[15,TRIALLowSalinity]"
[151] "r_FISHID[16,TRIALLowSalinity]"
[152] "r_FISHID[17,TRIALLowSalinity]"
[153] "r_FISHID[18,TRIALLowSalinity]"
[154] "r_FISHID[19,TRIALLowSalinity]"
[155] "r_FISHID[20,TRIALLowSalinity]"
[156] "r_FISHID[21,TRIALLowSalinity]"
[157] "r_FISHID[22,TRIALLowSalinity]"
[158] "r_FISHID[23,TRIALLowSalinity]"
[159] "r_FISHID[24,TRIALLowSalinity]"
[160] "r_FISHID[25,TRIALLowSalinity]"
[161] "r_FISHID[26,TRIALLowSalinity]"
[162] "r_FISHID[27,TRIALLowSalinity]"
[163] "r_FISHID[28,TRIALLowSalinity]"
[164] "r_FISHID[29,TRIALLowSalinity]"
[165] "r_FISHID[30,TRIALLowSalinity]"
[166] "r_FISHID[31,TRIALLowSalinity]"
[167] "r_FISHID[32,TRIALLowSalinity]"
[168] "r_FISHID[33,TRIALLowSalinity]"
[169] "r_FISHID[34,TRIALLowSalinity]"
[170] "r_FISHID[35,TRIALLowSalinity]"
[171] "r_FISHID[36,TRIALLowSalinity]"
[172] "r_FISHID[37,TRIALLowSalinity]"
[173] "r_FISHID[38,TRIALLowSalinity]"
[174] "r_FISHID[39,TRIALLowSalinity]"
[175] "r_FISHID[40,TRIALLowSalinity]"
[176] "r_FISHID[41,TRIALLowSalinity]"
[177] "r_FISHID[42,TRIALLowSalinity]"
[178] "r_FISHID[43,TRIALLowSalinity]"
[179] "r_FISHID[44,TRIALLowSalinity]"
[180] "r_FISHID[45,TRIALLowSalinity]"
[181] "r_FISHID[46,TRIALLowSalinity]"
[182] "r_FISHID[47,TRIALLowSalinity]"
[183] "r_FISHID[48,TRIALLowSalinity]"
[184] "r_FISHID[49,TRIALLowSalinity]"
[185] "r_FISHID[50,TRIALLowSalinity]"
[186] "r_FISHID[51,TRIALLowSalinity]"
[187] "r_FISHID[52,TRIALLowSalinity]"
[188] "r_FISHID[53,TRIALLowSalinity]"
[189] "r_FISHID[54,TRIALLowSalinity]"
[190] "r_FISHID[55,TRIALLowSalinity]"
[191] "r_FISHID[56,TRIALLowSalinity]"
[192] "r_FISHID[57,TRIALLowSalinity]"
[193] "r_FISHID[58,TRIALLowSalinity]"
[194] "r_FISHID[59,TRIALLowSalinity]"
[195] "r_FISHID[60,TRIALLowSalinity]"
[196] "prior_Intercept"
[197] "prior_b_TRIALHypoxia"
[198] "prior_b_TRIALLowSalinity"
[199] "prior_b_SMR_contr"
[200] "prior_b_MASS"
[201] "prior_sigma"
[202] "prior_sd_FISHID"
[203] "prior_cor_FISHID"
[204] "lprior"
[205] "lp__"
[206] "accept_stat__"
[207] "stepsize__"
[208] "treedepth__"
[209] "n_leapfrog__"
[210] "divergent__"
[211] "energy__"
# A tibble: 26,400 × 5
# Groups: .variable [11]
.chain .iteration .draw .variable .value
<int> <int> <int> <chr> <dbl>
1 1 1 1 b_Intercept 139.
2 1 2 2 b_Intercept 132.
3 1 3 3 b_Intercept 135.
4 1 4 4 b_Intercept 169.
5 1 5 5 b_Intercept 144.
6 1 6 6 b_Intercept 154.
7 1 7 7 b_Intercept 178.
8 1 8 8 b_Intercept 155.
9 1 9 9 b_Intercept 132.
10 1 10 10 b_Intercept 149.
# ℹ 26,390 more rows
We can then summarise this
# A tibble: 11 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_Intercept 146. 1.01e+2 189. 0.95 median hdci
2 b_MASS 0.118 -3.51e-1 0.573 0.95 median hdci
3 b_SMR_contr -19.7 -2.67e+1 -12.4 0.95 median hdci
4 b_TRIALHypoxia -108. -1.53e+2 -68.3 0.95 median hdci
5 b_TRIALHypoxia:SMR_contr 11.0 2.51e+0 18.6 0.95 median hdci
6 b_TRIALLowSalinity -72.9 -1.37e+2 -10.2 0.95 median hdci
7 b_TRIALLowSalinity:SMR_con… 5.89 -6.65e+0 17.7 0.95 median hdci
8 sd_FISHID__Intercept 11.8 6.08e+0 17.8 0.95 median hdci
9 sd_FISHID__TRIALHypoxia 6.07 2.69e-3 12.8 0.95 median hdci
10 sd_FISHID__TRIALLowSalinity 28.0 2.03e+1 35.5 0.95 median hdci
11 sigma 13.4 9.83e+0 17.4 0.95 median hdci
norin.brm4 |>
gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)Warning: `stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95), labels =
scales::percent_format()))` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()))` instead.
norin.brm4 |>
gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")norin.brm4 |>
gather_draws(`b_Intercept|b_.*`, regex = TRUE) |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()norin.brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")norin.brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
geom_vline(xintercept = 0, linetype = "dashed")norin.brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
geom_vline(xintercept = 0, linetype = "dashed")Picking joint bandwidth of 2.12
## Or in colour
norin.brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
geom_density_ridges_gradient(
aes(
x = exp(.value),
y = .variable,
fill = stat(x)
),
alpha = 0.4, colour = "white",
quantile_lines = TRUE,
quantiles = c(0.025, 0.975)
) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans()) +
scale_fill_viridis_c(option = "C")Picking joint bandwidth of 3.06
# A tibble: 2,400 × 214
.chain .iteration .draw b_Intercept b_TRIALHypoxia b_TRIALLowSalinity
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 139. -153. -73.1
2 1 2 2 132. -103. -73.4
3 1 3 3 135. -107. -39.0
4 1 4 4 169. -140. -81.1
5 1 5 5 144. -127. -3.04
6 1 6 6 154. -104. 1.02
7 1 7 7 178. -105. -47.8
8 1 8 8 155. -105. -55.8
9 1 9 9 132. -80.0 -61.9
10 1 10 10 149. -106. -109.
# ℹ 2,390 more rows
# ℹ 208 more variables: b_SMR_contr <dbl>, b_MASS <dbl>,
# `b_TRIALHypoxia:SMR_contr` <dbl>, `b_TRIALLowSalinity:SMR_contr` <dbl>,
# sd_FISHID__Intercept <dbl>, sd_FISHID__TRIALHypoxia <dbl>,
# sd_FISHID__TRIALLowSalinity <dbl>,
# cor_FISHID__Intercept__TRIALHypoxia <dbl>,
# cor_FISHID__Intercept__TRIALLowSalinity <dbl>, …
# A tibble: 2,400 × 75
.chain .iteration .draw b_Intercept b_TRIALHypoxia b_TRIALLowSalinity
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 139. -153. -73.1
2 1 2 2 132. -103. -73.4
3 1 3 3 135. -107. -39.0
4 1 4 4 169. -140. -81.1
5 1 5 5 144. -127. -3.04
6 1 6 6 154. -104. 1.02
7 1 7 7 178. -105. -47.8
8 1 8 8 155. -105. -55.8
9 1 9 9 132. -80.0 -61.9
10 1 10 10 149. -106. -109.
# ℹ 2,390 more rows
# ℹ 69 more variables: b_SMR_contr <dbl>, b_MASS <dbl>,
# `b_TRIALHypoxia:SMR_contr` <dbl>, `b_TRIALLowSalinity:SMR_contr` <dbl>,
# sd_FISHID__Intercept <dbl>, cor_FISHID__Intercept__TRIALHypoxia <dbl>,
# cor_FISHID__Intercept__TRIALLowSalinity <dbl>, Intercept <dbl>,
# `r_FISHID[1,Intercept]` <dbl>, `r_FISHID[2,Intercept]` <dbl>,
# `r_FISHID[3,Intercept]` <dbl>, `r_FISHID[4,Intercept]` <dbl>, …
Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
# A tibble: 2,400 × 205
b_Intercept b_TRIALHypoxia b_TRIALLowSalinity b_SMR_contr b_MASS
<dbl> <dbl> <dbl> <dbl> <dbl>
1 139. -153. -73.1 -19.9 0.339
2 132. -103. -73.4 -17.1 0.290
3 135. -107. -39.0 -18.5 0.336
4 169. -140. -81.1 -22.7 -0.00966
5 144. -127. -3.04 -19.4 0.135
6 154. -104. 1.02 -22.4 0.245
7 178. -105. -47.8 -23.9 -0.0865
8 155. -105. -55.8 -20.6 0.0615
9 132. -80.0 -61.9 -15.1 0.0266
10 149. -106. -109. -21.3 0.269
# ℹ 2,390 more rows
# ℹ 200 more variables: `b_TRIALHypoxia:SMR_contr` <dbl>,
# `b_TRIALLowSalinity:SMR_contr` <dbl>, sd_FISHID__Intercept <dbl>,
# sd_FISHID__TRIALHypoxia <dbl>, sd_FISHID__TRIALLowSalinity <dbl>,
# cor_FISHID__Intercept__TRIALHypoxia <dbl>,
# cor_FISHID__Intercept__TRIALLowSalinity <dbl>,
# cor_FISHID__TRIALHypoxia__TRIALLowSalinity <dbl>, sigma <dbl>, …
y ymin ymax .width .point .interval
1 0.4980791 0.4431186 0.5484889 0.95 median hdci
y ymin ymax .width .point .interval
1 0.6128526 0.5304293 0.6687602 0.95 median hdci
y ymin ymax .width .point .interval
1 0.8459696 0.7474859 0.9157748 0.95 median hdci
Region of Practical Equivalence
[1] 3.38454
# Proportion of samples inside the ROPE [-3.38, 3.38]:
Parameter | inside ROPE
----------------------------------------
Intercept | 0.00 %
TRIALHypoxia | 0.00 %
TRIALLowSalinity | 0.00 %
SMR_contr | 0.00 %
MASS | 100.00 %
TRIALHypoxia:SMR_contr | 0.70 %
TRIALLowSalinity:SMR_contr | 27.76 %
norin.brm4 |> modelsummary(
statistic = c("conf.low", "conf.high"),
shape = term ~ statistic,
exponentiate = FALSE
)Warning:
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
This warning appears once per session.
| (1) | |||
|---|---|---|---|
| Est. | 2.5 % | 97.5 % | |
| b_Intercept | 146.371 | 103.337 | 190.743 |
| b_TRIALHypoxia | -107.503 | -149.708 | -64.556 |
| b_TRIALLowSalinity | -72.901 | -134.697 | -5.921 |
| b_SMR_contr | -19.665 | -26.600 | -12.287 |
| b_MASS | 0.118 | -0.340 | 0.590 |
| b_TRIALHypoxia × SMR_contr | 10.981 | 2.857 | 18.973 |
| b_TRIALLowSalinity × SMR_contr | 5.895 | -6.628 | 17.750 |
| sd_FISHID__Intercept | 11.770 | 5.723 | 17.449 |
| sd_FISHID__TRIALHypoxia | 6.071 | 0.392 | 14.586 |
| sd_FISHID__TRIALLowSalinity | 27.978 | 20.799 | 36.427 |
| cor_FISHID__Intercept__TRIALHypoxia | -0.397 | -0.881 | 0.616 |
| cor_FISHID__Intercept__TRIALLowSalinity | 0.049 | -0.336 | 0.613 |
| cor_FISHID__TRIALHypoxia__TRIALLowSalinity | -0.550 | -0.952 | 0.519 |
| sigma | 13.382 | 9.501 | 17.208 |
| Num.Obs. | 180 | ||
| R2 | 0.846 | ||
| R2 Adj. | 0.750 | ||
| R2 Marg. | 0.498 | ||
| ELPD | -780.1 | ||
| ELPD s.e. | 9.2 | ||
| LOOIC | 1560.3 | ||
| LOOIC s.e. | 18.4 | ||
| WAIC | 1526.5 | ||
| RMSE | 9.26 | ||
| r2.adjusted.marginal | 0.308806833536584 | ||
9 Predictions / further analyses
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
contrast estimate lower.HPD upper.HPD
HighTemperature - Hypoxia -10.98 -18.94 -2.84
HighTemperature - LowSalinity -5.89 -17.83 6.37
Hypoxia - LowSalinity 4.85 -8.74 17.46
Point estimate displayed: median
HPD interval probability: 0.95
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
Conclusions:
- the partial slope associated with control SMR during High Temperature was found to be 10.98 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
- the partial slope associated with control SMR during High Temperature was found to be 5.89 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
- the partial slope associated with control SMR during Hypoxia was found to be 4.85 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.
$SMR_contr
Mean Lower Upper
5.178857 3.926671 6.431044
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
SMR_contr = 3.93:
contrast estimate lower.HPD upper.HPD
HighTemperature - Hypoxia 64.32 53.7 76.28
HighTemperature - LowSalinity 49.41 33.1 66.86
Hypoxia - LowSalinity -15.09 -33.4 4.29
SMR_contr = 5.18:
contrast estimate lower.HPD upper.HPD
HighTemperature - Hypoxia 50.90 45.0 56.15
HighTemperature - LowSalinity 41.71 33.3 50.35
Hypoxia - LowSalinity -8.89 -18.4 1.43
SMR_contr = 6.43:
contrast estimate lower.HPD upper.HPD
HighTemperature - Hypoxia 36.97 25.8 49.29
HighTemperature - LowSalinity 34.31 16.2 51.17
Hypoxia - LowSalinity -2.86 -23.2 15.96
Point estimate displayed: median
HPD interval probability: 0.95
norin.brm4 |>
emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
pairs() |>
gather_emmeans_draws() |>
median_hdci()Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 9 × 8
contrast SMR_contr .value .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 HighTemperature - Hypo… 3.93 64.3 53.7 76.3 0.95 median hdci
2 HighTemperature - Hypo… 5.18 50.9 45.0 56.1 0.95 median hdci
3 HighTemperature - Hypo… 6.43 37.0 25.8 49.3 0.95 median hdci
4 HighTemperature - LowS… 3.93 49.4 33.1 66.9 0.95 median hdci
5 HighTemperature - LowS… 5.18 41.7 33.3 50.4 0.95 median hdci
6 HighTemperature - LowS… 6.43 34.3 16.5 51.7 0.95 median hdci
7 Hypoxia - LowSalinity 3.93 -15.1 -33.6 4.15 0.95 median hdci
8 Hypoxia - LowSalinity 5.18 -8.89 -18.6 1.33 0.95 median hdci
9 Hypoxia - LowSalinity 6.43 -2.86 -23.2 16.0 0.95 median hdci
norin.brm4 |>
emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
pairs() |>
tidy_draws() |>
summarise_draws(median)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 9 × 2
variable median
<chr> <dbl>
1 contrast HighTemperature - Hypoxia SMR_contr 5.17885746916438 50.9
2 contrast HighTemperature - LowSalinity SMR_contr 5.17885746916438 41.7
3 contrast Hypoxia - LowSalinity SMR_contr 5.17885746916438 -8.89
4 contrast HighTemperature - Hypoxia SMR_contr 3.92667075536741 64.3
5 contrast HighTemperature - LowSalinity SMR_contr 3.92667075536741 49.4
6 contrast Hypoxia - LowSalinity SMR_contr 3.92667075536741 -15.1
7 contrast HighTemperature - Hypoxia SMR_contr 6.43104418296134 37.0
8 contrast HighTemperature - LowSalinity SMR_contr 6.43104418296134 34.3
9 contrast Hypoxia - LowSalinity SMR_contr 6.43104418296134 -2.86
Conclusions:
- as control SMR increases, the magnitude of differences between the trials becomes reduced.
- the High Temperature trial resulted in the greatest positive change in SMR.
- the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.
Of course, as with other analyses, it is possible to extract the entire posterior and use it to calculate exceedence probabilities etc.
norin.em <- norin.brm4 |>
emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
pairs() |>
gather_emmeans_draws() |>
mutate(Fit = .value)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
# A tibble: 21,600 × 7
# Groups: contrast, SMR_contr [9]
contrast SMR_contr .chain .iteration .draw .value Fit
<fct> <dbl> <int> <int> <int> <dbl> <dbl>
1 HighTemperature - Hypoxia 5.18 NA NA 1 49.6 49.6
2 HighTemperature - Hypoxia 5.18 NA NA 2 52.5 52.5
3 HighTemperature - Hypoxia 5.18 NA NA 3 55.3 55.3
4 HighTemperature - Hypoxia 5.18 NA NA 4 52.4 52.4
5 HighTemperature - Hypoxia 5.18 NA NA 5 49.7 49.7
6 HighTemperature - Hypoxia 5.18 NA NA 6 50.4 50.4
7 HighTemperature - Hypoxia 5.18 NA NA 7 51.3 51.3
8 HighTemperature - Hypoxia 5.18 NA NA 8 50.0 50.0
9 HighTemperature - Hypoxia 5.18 NA NA 9 55.3 55.3
10 HighTemperature - Hypoxia 5.18 NA NA 10 52.6 52.6
# ℹ 21,590 more rows
# A tibble: 3 × 7
contrast Fit .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 HighTemperature - Hypoxia 50.9 28.9 73.2 0.95 median hdci
2 HighTemperature - LowSalinity 41.8 22.2 62.6 0.95 median hdci
3 Hypoxia - LowSalinity -8.92 -30.0 11.6 0.95 median hdci
# A tibble: 9 × 8
contrast SMR_contr Fit .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 HighTemperature - Hypo… 3.93 64.3 53.7 76.3 0.95 median hdci
2 HighTemperature - Hypo… 5.18 50.9 45.0 56.1 0.95 median hdci
3 HighTemperature - Hypo… 6.43 37.0 25.8 49.3 0.95 median hdci
4 HighTemperature - LowS… 3.93 49.4 33.1 66.9 0.95 median hdci
5 HighTemperature - LowS… 5.18 41.7 33.3 50.4 0.95 median hdci
6 HighTemperature - LowS… 6.43 34.3 16.5 51.7 0.95 median hdci
7 Hypoxia - LowSalinity 3.93 -15.1 -33.6 4.15 0.95 median hdci
8 Hypoxia - LowSalinity 5.18 -8.89 -18.6 1.33 0.95 median hdci
9 Hypoxia - LowSalinity 6.43 -2.86 -23.2 16.0 0.95 median hdci
## norin.em |>
## group_by(contrast) |>
## summarize(P=sum(Fit>0)/n())
norin.em |>
group_by(contrast, SMR_contr) |>
summarise(
P = mean(Fit > 0),
P2 = 1 - P
)`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.
# A tibble: 9 × 4
# Groups: contrast [3]
contrast SMR_contr P P2
<fct> <dbl> <dbl> <dbl>
1 HighTemperature - Hypoxia 3.93 1 0
2 HighTemperature - Hypoxia 5.18 1 0
3 HighTemperature - Hypoxia 6.43 1 0
4 HighTemperature - LowSalinity 3.93 1 0
5 HighTemperature - LowSalinity 5.18 1 0
6 HighTemperature - LowSalinity 6.43 1.000 0.000417
7 Hypoxia - LowSalinity 3.93 0.0542 0.946
8 Hypoxia - LowSalinity 5.18 0.0342 0.966
9 Hypoxia - LowSalinity 6.43 0.397 0.603
10 Summary figures
norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.brm4 |>
emmeans(~ SMR_contr | TRIAL, at = norin.grid) |>
as.data.frame()Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
SMR_contr TRIAL emmean lower.HPD upper.HPD
3.984319 HighTemperature 73.88891 63.57026 83.21900
4.011756 HighTemperature 73.36914 63.20446 82.49482
4.039192 HighTemperature 72.84601 62.88133 81.85916
4.066629 HighTemperature 72.28818 63.18352 81.79213
4.094065 HighTemperature 71.74851 62.75470 80.98760
4.121502 HighTemperature 71.22935 62.35300 80.23893
Point estimate displayed: median
HPD interval probability: 0.95
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
theme_classic() +
theme(
legend.position = c(0.99, 0.99),
legend.justification = c(1, 1)
)Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## The .fixed values are the predicted values without random effects
obs <- norin.brm4 |>
augment() |>
mutate(PartialObs = .fitted + .resid)
ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
## geom_point(data=norin, aes(y=CHANGE), color='gray') +
theme_classic()g1 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = obs, aes(y = PartialObs, color = TRIAL)) +
theme_classic()
g2 <- ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
geom_line(aes(, color = TRIAL)) +
geom_point(data = norin, aes(y = CHANGE, color = TRIAL)) +
theme_classic()
g1 + g211 References
norin <- norin |> mutate(
FISHID = factor(FISHID),
TRIAL = factor(TRIAL)
)
ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth(method = "lm") +
geom_point()
# ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
# geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
geom_point() +
geom_line()
# ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = MASS, color = TRIAL)) +
geom_point() +
geom_smooth(method = "lm")
norin.rstanarm <- stan_glmer(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
data = norin,
prior_PD = TRUE,
iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0
)
prior_summary(norin.rstanarm)
posterior_vs_prior(norin.rstanarm,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y"), pars = c("(Intercept)")
)
ggpredict(norin.rstanarm, ~ TRIAL * SMR_contr) |> plot(show_data = TRUE)
norin.rstanarm |> get_variables()
plot(norin.rstanarm, "mcmc_trace", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_acf_bar", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_rhat_hist", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
plot(norin.rstanarm, "mcmc_neff_hist", regex_pars = "^.Intercept|TRIAL|SMR|MASS|[sS]igma")
# norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
# iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 <- stan_glmer(CHANGE ~ (SMR_contr | FISHID) + TRIAL * SMR_contr + MASS,
data = norin,
prior_PD = FALSE,
iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0, cores = 3
)
norin.rstanarm1 <- update(norin.rstanarm1, prior_PD = FALSE)
norin.rstanarm2 <- stan_glmer(CHANGE ~ (TRIAL * SMR_contr | FISHID) + TRIAL * SMR_contr + MASS,
data = norin,
prior_PD = FALSE,
iter = 5000, warmup = 2000, chains = 3, thin = 5, refresh = 0, cores = 3
)
posterior_vs_prior(norin.rstanarm1,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y"), pars = c("(Intercept)")
)
ggpredict(norin.rstanarm1, ~ TRIAL * SMR_contr) |> plot(show_data = TRUE)
norin.rstanarm1 |> get_variables()
plot(norin.rstanarm1, "mcmc_trace", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_acf_bar", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_rhat_hist", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
plot(norin.rstanarm1, "mcmc_neff_hist", regex_pars = "^.Intercept|TRIAL|^SMR|MASS|[sS]igma")
(l.1 <- loo(norin.rstanarm))
(l.2 <- loo(norin.rstanarm1))
loo_compare(l.1, l.2)
preds <- posterior_predict(norin.rstanarm, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median)
)
plot(norin.resids)
g <- ggpredict(norin.rstanarm) |> plot()
do.call("grid.arrange", g)
# ggemmeans(norin.rstan, ~TRIAL)
summary(norin.rstanarm)
nms <- norin.rstanarm1 |> get_variables()
wch <- grep("^.Intercept|TRIAL|^SMR|[sS]igma", nms)
tidyMCMC(norin.rstanarm$stanfit,
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE, pars = nms[wch], estimate.method = "median"
)
tidyMCMC(norin.rstanarm1$stanfit,
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE, pars = nms[wch], estimate.method = "median"
)
norin.grid <- with(norin, list(SMR_contr = seq(min(SMR_contr), max(SMR_contr), len = 100)))
newdata <- emmeans(norin.rstanarm, ~ TRIAL | SMR_contr, at = norin.grid) |> as.data.frame()
head(newdata)
ggplot(newdata, aes(y = emmean, x = SMR_contr, color = TRIAL)) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3, color = NA) +
geom_line()
norin.grid <- with(norin, list(SMR_contr = c(min(SMR_contr), mean(SMR_contr), max(SMR_contr))))
emmeans(norin.rstan, pairwise ~ TRIAL | SMR_contr, at = norin.grid)
norin.em <- emmeans(norin.rstan, pairwise ~ TRIAL | SMR_contr, at = norin.grid)$contrast |>
gather_emmeans_draws() |>
mutate(Fit = .value)
norin.em
norin.em |>
group_by(contrast) |>
median_hdci(Fit)
norin.em |>
group_by(contrast, SMR_contr) |>
median_hdci(Fit)
## norin.em |>
## group_by(contrast) |>
## summarize(P=sum(Fit>0)/n())
norin.em |>
group_by(contrast, SMR_contr) |>
summarize(P = mean(Fit > 0))
bayes_R2(norin.rstanarm, re.form = NA) |> median_hdi()
bayes_R2(norin.rstanarm, re.form = ~ (1 | FISHID)) |> median_hdi()
# bayes_R2(norin.rstan1, re.form=~(SMR_contr|FISHID)) |> median_hdi11.1 brms
norin <- norin |> mutate(
FISHID = factor(FISHID),
TRIAL = factor(TRIAL)
)
ggplot(norin, aes(y = CHANGE, x = TRIAL)) +
geom_boxplot()
ggplot(norin, aes(y = CHANGE, x = SMR_contr, shape = TRIAL, color = TRIAL)) +
geom_smooth(method = "lm") +
geom_point()
ggplot(norin, aes(y = CHANGE, x = MASS, shape = TRIAL, color = TRIAL)) +
geom_smooth(method = "lm") +
geom_point()
ggplot(norin, aes(y = CHANGE, x = as.numeric(FISHID), color = TRIAL)) +
geom_point() +
geom_line()
## ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
## ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')
norin |>
group_by(TRIAL) |>
summarise(
median(CHANGE),
mad(CHANGE)
)
priors <- prior(normal(50, 20), class = "Intercept") +
prior(normal(0, 10), class = "b") +
prior(gamma(2, 1), class = "sigma") +
prior(gamma(2, 1), class = "sd")
norin.form <- bf(CHANGE ~ (1 | FISHID) + TRIAL * SMR_contr + MASS,
family = gaussian
)
norin.brm1 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = "yes",
iter = 5000, warmup = 2000,
chains = 3, thin = 5, refresh = 0
)
norin.brm1 |> get_variables()
pars <- norin.brm1 |> get_variables()
wch <- grepl("^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd", pars, perl = TRUE)
stan_trace(norin.brm1$fit, pars = pars[wch])
stan_ac(norin.brm1$fit, pars = pars[wch])
stan_rhat(norin.brm1$fit, pars = pars[wch])
stan_ess(norin.brm1$fit, pars = pars[wch])
## mcmc_plot(norin.brms, type='trace',
## regex_pars='^b.*|sigma|^sd')
## mcmc_plot(norin.brms, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms, type='acf_bar',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms, type='rhat_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brms, type='neff_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
preds <- posterior_predict(norin.brm1, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median),
integerResponse = FALSE
)
plot(norin.resids)
# norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
# iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.form <- bf(CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS,
family = gaussian
)
norin.brm2 <- brm(norin.form,
data = norin,
prior = priors,
sample_prior = "yes",
iter = 5000, warmup = 2000,
chains = 3, thin = 5, refresh = 0, cores = 3,
control = list(adapt_delta = 0.99)
)
norin.brm2 |> get_variables()
pars <- norin.brm2 |> get_variables()
## wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)
wch <- grepl("^b_.*|[sS]igma|^sd_.*", pars, perl = TRUE)
stan_trace(norin.brm2$fit, pars = pars[wch])
stan_ac(norin.brm2$fit, pars = pars[wch])
stan_rhat(norin.brm2$fit) # , pars=pars[wch])
stan_ess(norin.brm2$fit) # , pars=pars[wch])
## mcmc_plot(norin.brm2, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2, type='trace',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2, type='acf_bar',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2, type='rhat_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
## mcmc_plot(norin.brm2, type='neff_hist',
## regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
(l.1 <- loo(norin.brm1))
(l.2 <- loo(norin.brm2))
loo_compare(l.1, l.2)
preds <- posterior_predict(norin.brm2, nsamples = 250, summary = FALSE)
norin.resids <- createDHARMa(
simulatedResponse = t(preds),
observedResponse = norin$CHANGE,
fittedPredictedResponse = apply(preds, 2, median)
)
plot(norin.resids)
g <- norin.brm2 |>
conditional_effects() |>
plot(points = TRUE, ask = FALSE)
library(patchwork)
g[[1]] + g[[2]] + g[[3]] + g[[4]]
## g=ggpredict(norin.brms1) |> plot
## library(patchwork)
## g[[1]] + g[[2]] + g[[3]]
## do.call('grid.arrange', g)
ggemmeans(norin.brm2, ~TRIAL) |> plot()
summary(norin.brm2)
tidyMCMC(norin.brm2$fit,
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE, estimate.method = "median"
) |>
slice(1:11)
pars <- norin.brm2 |> get_variables()
wch <- grep("^b.Intercept|TRIAL|^b.*SMR|[sS]igma|^sd", pars)
tidyMCMC(norin.brms1$fit,
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE, pars = pars[wch], estimate.method = "median"
)
bayes_R2(norin.brm2, re.form = NA, summary = FALSE) |>
median_hdci()
bayes_R2(norin.brm2, re.form = ~ (1 | FISHID), summary = FALSE) |>
median_hdci()
bayes_R2(norin.brm2, re.form = ~ (TRIAL | FISHID), summary = FALSE) |>
median_hdci()
emmeans(norin.brm2, pairwise ~ TRIAL)
norin.em <- norin.brm2 |>
emmeans(~TRIAL) |>
pairs() |>
gather_emmeans_draws() |>
mutate(Fit = .value)
norin.em |>
group_by(contrast) |>
median_hdi()
norin.em |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = contrast,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE) +
theme_bw()
norin.em |>
group_by(contrast) |>
summarize(P = mean(Fit > 0))
norin.grid <- with(norin, list(SMR_contr = c(
min(SMR_contr),
mean(SMR_contr),
max(SMR_contr)
)))
norin.em <- norin.brm2 |>
emmeans(~ TRIAL | SMR_contr, at = norin.grid) |>
pairs() |>
gather_emmeans_draws()
norin.em |> head()
norin.em |>
group_by(contrast, SMR_contr) |>
median_hdi()
norin.em |>
group_by(contrast, SMR_contr) |>
summarize(P = mean(.value > 0))
norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.brm2 |>
emmeans(~ SMR_contr | TRIAL, at = norin.grid) |>
as.data.frame()
head(newdata)
partial.obs <- norin |>
mutate(
Pred = predict(norin.brm2, re.form = NA, summary = TRUE)[, "Estimate"],
Resid = resid(norin.brm2)[, "Estimate"],
Obs = Pred + Resid
)
ggplot(newdata, aes(y = emmean, x = SMR_contr, color = TRIAL)) +
geom_point(data = partial.obs, aes(y = Obs)) +
## geom_point(data=partial.obs, aes(y=CHANGE), shape=2) +
geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3, color = NA) +
geom_line()